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Abstract 



The neutralino, the lightest stable supersymmetric particle, is a strong theoretical candi- 
date for the missing astronomical "dark matter". Depending on their annihilation cross 
section, relic neutralinos from early formation of the Universe trapped in orbits around 
massive objects may currently be annihilating at measurable rates. The Minimal Su- 
persymmetric extension of the Standard Model predicts that the gamma rays emerging 
from one of the annihilation modes will give a distinct monochromatic signal with en- 
ergy between lOOGeV and lOTeV, depending on the neutralino mass. An additional 
"continuum" spectrum signal of photons will be produced by the decay of secondaries 
produced in the non-photonic annihilation modes. 

Milagro is an air shower array which uses the water Cherenkov technique and is 
capable of detecting TeV gamma rays from the direction of the Sun with an angular 
resolution of less than a degree. It is the first instrument capable of establishing a limit 
on the gamma-ray flux from neutralino annihilations near the Sun. 

In this report results of a search for neutralino to photon annihilation with the Mi- 
lagro gamma-ray observatory are presented. Results of a Monte Carlo computer sim- 
ulation of the neutralino annihilation density in the Solar System suggest that a large 
portion of neutralino annihilations (40 — 50%) happens outside the Sun which may give 
rise to a detectable gamma-ray signal from the solar region. No significant gamma-ray 
signal was observed from the Sun resulting in an upper limit on the sought for photon 
flux. The upper limit can be translated to a neutralino-mass dependent limit on the prod- 
uct of the neutralino -proton scattering crossection a px , the integrated photon yield per 
neutralino in neutralino-neutralino annihilation 6 7 and the local galactic halo dark mat- 
ter density p . For example, assuming a ITeV neutralino and ignoring the continuum 
contribution to the signal gives an upper limit of a3 {G eVcm-*) m=^k^) K < 2 - 3 - 
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Chapter 1 
Introduction 



Chapter 1, like Chapter 3, had not actually 
told him anything that he did not know; it 
had merely systematized the knowledge that 
he possessed already. 

George Orwell, "1984" 



1.1 The dark matter problem. 

Perhaps, there is no problem of greater importance to cosmology and astrophysics than 
that of the "dark matter". It is centered around the notion that there may exist an enor- 
mous amount of non-luminous matter in the Universe. The presence of the matter, 
which does not radiate and can not be seen directly, can only be inferred by observing 
the effects it has on other directly observable astronomical objects. 

It has always been known that there is matter in the sky which does not emit any 
kind of radiation. For instance, the planets do not shine, but their contribution to the 
mass of the solar system is negligible, so worrying about non-luminous matter was not 
of a great concern. 

The first evidence that there is a significant amount of dark matter came from Zwicky 
[53], in the thirties, from investigations of clusters of galaxies. It was found that veloci- 
ties of the galaxies in a cluster were about 10 times larger than expected, indicating that 
there is invisible gravitating matter in a cluster, holding the galaxies together. 

Somewhat more reliable evidence was found in the 1970s by Rubin [44] by studying 
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Figure 1.1: A "typical" rotation curve of a "typical" galaxy, determined from 21cm 
observations. [45] 

the rotation curves of spiral galaxies. Kepler's law states that the rotational velocity 
around a gravitational center depends only on the distance from the center and on the 
total mass contained within the orbit. Thus, one expects: 

GM(r) = v 2 r 

where v - is the circular velocity at a distance r from the center of the galaxy, M(r) is 
the mass enclosed in the sphere r and spherical symmetry is assumed. 

If the mass were associated with light (luminous matter) v would decrease as r -1 / 2 
beyond the point where the light cuts off. However, it was found (See, for instance [45, 
41], or fig 1.1) that v ~ const, corresponding to M(r) oc r, implying existence of dark 
halos around spiral galaxies. The halos could be made of brown dwarf stars, jupiters, 
planets and 1OOM black holes. Collectively, these objects are called MACHOs 1 and 
are the main baryonic dark matter candidates. 

Dark matter has important consequences for the evolution of the Universe. The stan- 
dard, Hot Big Bang cosmology is remarkably successful: it provides a reliable and tested 

'MaCHO is Massive Compact Halo Object 



2 



account of the history of the Universe from at least t ~ 10~ 2 s until today (t ~ 14 Gyr). 
At present, there is no strong experimental evidence contradicting the theory. According 
to the theory, the Universe must conform to one of three possible types with negative, 
positive or zero curvature. The value of the cosmological density parameter, 2 fltotai, 
determines which of the three possibilities applies to our world. There is, however, a 
somewhat philosophical or even aesthetical argument that makes VL to tai = 1 attractive. 
The point is that as the Universe evolves, the value of VL total changes. In fact, the value 
of Vt total = 1 is unstable. If the Universe is open VL total < 1, it will expand forever, until 
it is totally empty VL 0. On contrary, if it is closed n to tai > 1, it will recollapse to 
a state with extremely high density f2 oo. The inflationary cosmology [31], which 
provides the most compelling explanation for the smoothness of the Cosmic Microwave 
Background Radiation (CMBR), predicts that the early Universe was extremely close to 
flat \Vltotai — 1| < 0(1CT 60 ), leading to the belief that Vl tota i is exactly one. 

In fact, the most recent results from the studies of the CMBR with the WMAP [7] 
observatory yield VL tota i = 1-02 ± 0.02. The same study implies that the matter compo- 
nent of the total energy density is Vl m = 0.27 ±0.04 while ordinary baryonic component 
constitutes only about 15% of all matter in the Universe VL b = 0.044 ± 0.004. The rest 
of the energy density Vt\ = (Vt tota i — Vt b ) = 0.73 ± 0.04 is an unknown from of energy 
(so-called "dark energy"). 

In any event, the abundance of baryons is not likely to account for all matter even 
if fltotai turns out to be slightly less than unity and non-baryonic dark matter is almost 
required to dominate the Universe. The particles or fields which comprise nonbaryonic 
dark matter must have survived from the Big Bang, and therefore, must either be stable 
or have lifetimes in excess of the current age of the Universe. Among the non-baryonic 
dark matter candidates there are massive neutrinos, axions [49] and stable supersym- 
metric particles. 



2 £1 = p/pc, p - energy density of the Universe, p c - critical parameter; 

{< 1 negative curvature, the Universe will expand forever 
= 1 zero curvature Universe 
> 1 positive curvature, the Universe will recollapse, eventually 
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1.2 Supersymmetry. 



The main goal of the elementary particle physics is to devise a model which combines 
all particles and their interactions into one theory. The hope is that the development of 
supersymmetric theories (See, for example, [47]) is a step towards the stated goal. In 
these theories, bosonic and fermionic fields are allowed to transform into one another, 
and each particle is described by a multiplet containing bosons and fermions. In such 
models, loops, divergent in quantum field theories, cancel. Theoretical strong points of 
supersymmetry have motivated many accelerator searches for supersymmetric particles. 
Most of these have been guided by the Minimal Supersymmetric extension of the Stan- 
dard Model (MSSM) and are based on a missing-energy signature caused by the escape 
of the lightest supersymmetric particles. In the MSSM, the convergence of the renor- 
malized gauge couplings at the grand unification scale requires all masses of supersym- 
metric particles to appear between 100 (GeV) and 10 (TeV) [2]. Laboratory searches 
have set lower mass limits, requiring lightest supersymmetric particles in MSSM to pos- 
sess masses greater than 20 — 30 (GeV) [11]. Even though no convincing evidence for 
existence of supersymmetric particles has been found, they all have been given names. 
Bosonic ordinary particles have fermionic superpartners with the same name except 
with the suffix "ino" added, while fermionic ordinary particles have bosonic superpart- 
ners with prefix "s" added. For example, Higgsino is a superpartner for Higgs boson 
and selectron is a superpartner for electron. There are several superpartners which have 
the same quantum numbers and so can mix together in linear combinations. Since those 
do not necessarily correspond to any ordinary particle, they are given different names. 
For instance, the photino, Higgsino and Zino can mix into arbitrary combinations called 
the neutralinos. The lightest neutralino is a stable supersymmetric particle and makes 
the "best" candidate for a solution of the "dark matter problem"(first suggested in [42], 
also see [29] for an extensive review). 



1.3 Detection Methods. 

There are several ways to test the hypothesis that stable neutralinos exist and contribute 
to the dark matter. They include direct searches with extremely sensitive devices which 
can detect energy deposited by an elastically scattered neutralino off a nucleus and indi- 
rect searches which look for products of neutralino-neutralino annihilations. 
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WIMP Mass (GeV/c 2 ) 

Figure 1.2: Exclusion regions for neutralino-nucleon crossections obtained from differ- 
ent direct- search experiments. Closed contour is allowed region at 3<r confidence level 
from the DAMA experiment. The plot is adopted from [8]. 

1.3.1 Direct detection. 

The Italian/Chinese collaboration (DAMA) has reported an annual modulation in the to- 
tal count rate over 4 years. They interpret this as consistent with the annual modulation 
predicted for WIMPs [10]. This, however, is not a widely accepted result because of 
some possible modulating systematic errors. The CDMS experiment has obtained data 
that appear to exclude the DAMA result [1]. They reach a spin-independent WIMP- 
nucleon scrossection limit around 2 ■ 1CT 42 (cm 2 ) in the mass range 20 — 100 (GeV). 
Edelweiss has also released results that significantly cut into the DAMA allowed re- 
gion [8]. The summary of the limits of direct searches is shown on the figure 1.2. 
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1.3.2 Indirect Detection. 



Indirect searches also have received considerable attention from experimenters. For 
example, the Kamiokande and SuperKamiokande underground neutrino detectors have 
set limits on solar and terrestrial neutralino-induced muon fluxes [37, 25]. 

Another possible method for detecting dark matter particles is from their annihilation 
into 7-rays. One of the many possible gamma-producing channels is production of 
monochromatic gamma-rays: 

XX -> 77, XX -> %1 

Even though it is difficult to estimate the rates of these processes because of uncer- 
tainties in the supersymmetric parameters, cross sections and the neutralino distribution, 
since the annihilating neutralinos move at galactic velocities v/c ~ 1CT 3 the outgoing 
photons will give very distinct monochromatic signals 3 in each annihilation mode: 




which has no conceivable origin from any known astrophysical sources. 

As was mentioned earlier, neutralinos, if they are to be the dark matter, should have 
non-zero relic abundance today, but their number density is so small that almost no 
annihilations happen. An observation of such an event from some random point in the 
Universe is not feasible. However, since the density of neutralinos in the vicinity of a 
gravitational center will be larger than in other parts of the Universe and because the 
annihilation rate is proportional to the square of the neutralino density, there will be an 
enhanced flux of high energy 7-rays from such regions. Therefore, it is tempting to look 
at signals from well studied gravitating objects, such as nearby galaxies and the Milky 
Way Galaxy, and examine the energy spectrum for a monoenergetic signal. 

The present high energy gamma-ray experiments, such as EGRET and the Whip- 
ple atmospheric Cherenkov Telescope, lack the sensitivity to detect annihilation line 
fluxes predicted for most of the allowed supersymmetric models and halo profiles. How- 
ever, the next generation ground-based and satellite gamma-ray experiments, such as 

3 If these two lines can be resolved, the relative strength of the two could give a handle on the compo- 
sition of the neutralino. This is because despite the fact that the two processes are closely related, there 
are some differences which depend on the composition. 
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GRANITE-III, VERITAS and GLAST, will allow exploration of large portions of the 
MSSM parameter space, assuming that the dark matter density is peaked at the galactic 
center. [9] 

The Sun is also a large gravitating object and one could study the solar spectrum 
for the neutralino annihilation signal. Of course, that is possible only with a non- 
optical high resolution instrument, capable of monitoring the Sun at energies between 
100 (GeV) and 10 (TeV). Several semi-analytical estimates for the detection rates for 
several ground-based and satellite experiments are available in the literature. However, 
a more careful computer simulation following the decaying 3-D neutralino orbits with 
detailed elastic scattering and planetary perturbations accompanied by simulations of 
the solar magnetic field smearing and shadowing of the galactic cosmic rays by the Sun 
will provide a more definitive prediction on the neutralino annihilation rate. 

The structure of this work is the following: chapter 2 discusses how high energy 
cosmic particles can be detected. This is followed by a brief description of the Milagro 
detector, capable of monitoring the overhead sky at energies near 1 {TeV), in chapter 3. 
A presentation of the data analysis techniques employed in the current work is given in 
chapter 4. Chapter 5 discusses the computer simulations which are used to predict the 
gamma ray flux from the near-solar neutralino annihilations. Chapter 6 discusses the 
results of the search for the relic neutralinos and is followed by a summary in chapter 7. 



7 



Chapter 2 

Extensive Air Showers 



All one knew was that every quarter astro- 
nomical numbers of boots were produced on 
paper, while perhaps half the population of 
Oceania went barefoot. 

George Orwell "1984" 

There are several main reasons which govern the choice of a detector type to be used 
in high energy photon search from the Sun. First of all, it should be a non-optical de- 
vice capable of monitoring the solar region. Because, the Earth's atmosphere is opaque 
to gamma rays, satellite-based detectors need to be constructed to detect gamma rays. 
Indeed, small detectors sensitive to gamma rays at energies below a few GeV have been 
constructed and used successfully. 1 These detectors employ techniques developed for 
accelerator experiments where an incoming photon's direction is determined by e + e~ 
tracking detectors and the photon's energy is usually measured by a total-absorption 
calorimeter. However, the expected low and rapidly decreasing with photon energy 7- 
ray flux requires detectors with rather large collection areas and long exposure periods. 
Such detectors can be built on the surface of the Earth only. 

Even though direct detection of 7-rays is not possible by ground-based instruments, 
at energies above several GeV indirect gamma-ray detection is possible. Such very 
high energy photons initiate extensive air shower (EAS) cascades of secondary particles 
which are detectable by ground-based detectors. Knowledge of the EAS structure is 



'Future satellite detectors such as GLAST should register particles with energies as high as 
300 (GeV) [21]. 
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required to infer information about the primary photon. 

2.1 Development of EAS. 

Although an extrapolation from known particle physics might be necessary to describe 
the initial phase of the shower development, it is believed that the structure of the EAS is 
well understood. A high-energy primary photon interacts with electromagnetic fields of 
air molecules in the upper atmosphere producing an electron-positron pair which in turn 
produces high-energy photons via bremsstrahlung. The resulting electro-magnetic cas- 
cade grows geometrically as it propagates through the atmosphere. The shower growth 
stops when the mean energy of electrons and positrons falls below the critical energy 
(E c ~ 85 (MeV)) where the ionization energy-loss mechanism, which does not pro- 
duce additional shower particles, becomes dominant. After this point (called the shower 
maximum) the energy of particles and their number in the shower start to decrease as 
the shower continues its propagation towards the ground level. Nevertheless, a large 
number of shower particles may reach the ground and may be detected. 

Moreover, because the secondary particles are ultra-relativistic, they retain the di- 
rectionality of the incident gamma ray and the cascade arrives to the ground as a thin 
front perpendicular to the direction of the primary photon. The density of shower parti- 
cles in the front will decrease with distance from the extrapolated incident gamma-ray 
trajectory. This trajectory is called the core of the shower. 

The shower development is a stochastic process and while some analytical calcula- 
tions have been performed, computer simulations are generally employed to study the 
properties of the air shower cascades. 

2.1.1 Longitudinal Development of Extensive Air Showers. 

The average number of electrons N e and photons iV 7 in an electromagnetic shower 
can only depend on the primary energy E and the thickness of the traversed matter 
t. Moreover, if Eq is expressed in units of critical energy E c and t is in units of radiation 
lengths X , the number of electrons and photons is almost independent of the specific 
shower propagation medium. Usually, however, detectors can register particles with en- 
ergies above some E th , thus, often, it is desired to know the number of particles in a 
shower with energies greater than E th . According to [46] the average number of parti- 
cles N k (E , E th , t) of type k with energy above E th at atmospheric depth t in a shower 
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Eth, 


k 


=electron 


k=photon 


MeV 


A 


a 


b 


A 


a 


b 


1 


0.92 


0.00 


0.45 


4.80 


-0.88 


0.83 


5 


0.75 


0.19 


-1.22 


2.98 


-0.69 


-1.49 


10 


0.63 


0.35 


-2.57 


2.13 


-0.57 


-3.45 


20 


0.50 


0.53 


-4.22 


1.45 


-0.36 


-5.51 



Table 2. 1 : Values of parameters A, a and b for modified Greisen and NKG formulae. 
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Figure 2.1: Longitudinal development of electron (left) and photon (right) compontents 
of gamma-ray showers with E th = 1 (MeV) particle detection threshold. 



10 




Figure 2.2: Density of electrons (left) and photons (right) in a gamma-ray shower at 
atmospheric depth of 20 radiation lengths as a function of the core distance with E th = 
1 (MeV). Curves are normalized to the total number of respective particles. 



initiated by a photon with energy E can be described by a modified Greisen formula: 

N k (E ,E th ,t) = AkiE^^e^ 1 - 1 - 5 ^ 

y = ln^, t k = t + a k (E t h), s k 



E c ' ^ uin t k + 2y 

The parameterization is valid for 4 < t < 24 and 0.1 < E < 10 3 (TeV). The 
coefficients A k (E t h) and a k (E t h) are given in the table 2.1 The graphical illustration of 
the number of particles in a shower is presented in figure 2.1. 



2.1.2 Lateral Development of Extensive Air Showers. 

The average surface density p k (E , E th) t, r) of particles of type k in the shower front 
with energies greater than E th at a distance r from the shower axis and at the atmospheric 
depth t can be described by a modified Nishimura-Kamata-Greisen (NKG) function [46] 

/ 771 771 \ N k(E ,E th ,t) 

p k (E ,E th ,t,r) = -2 f(r/R k ,s k ) 
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Figure 2.3: Average arrival time of the shower front as a function of core distance (il- 
lustration). 



k t + b k (E th ) + 2y> I[X ' Z) 2vr B(zA-5-zf [i + X) 

where B(z, w) is the beta-function so that 2tt J °° f(x, z)xdx = 1. The characteristic 

scattering length for photons _R 7 = meC ^ 4n ^ a -X Q is the Moliere scattering unit and for 
electrons — R e = _R 7 /2. 

The values of parameters b(E th ) are given in the table 2.1 and the average density of 
photons and electrons per unit area is shown of figure 2.2 as a function of distance from 
the shower axis. 



2.1.3 Temporal Distribution of Extensive Air Shower Particles. 

Because the air shower detectors determine the primary particle direction using particle 
arrival times, knowledge of the shape of the shower front is important for achieving the 
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best possible angular resolution. Results of Monte Carlo simulations [46] of shower 
front are shown on figure 2.3. The shower front appears to have a parabolic shape as 
a function of distance from the core. At large atmospheric depths air shower photons 
travel faster than air shower electrons thus the photonic front is curved less than the 
electronic one. The thickness of the shower is defined by the distribution of the shower 
particle arrival times at distance r from the shower core and increases with the core 
distance r. At small core distances the fluctuations of arrival time around its average ap- 
pear to be smaller for the photon component than for the electron one and, consequently, 
photonic thickness is smaller than the electronic one. At large core distances, however, 
the electronic contribution to the shower is quite small due to electron ionization losses 
in the atmosphere compared to the photonic component. 



2.2 Cosmic rays. 



Among the particles which enter the Earth's atmosphere gamma rays present a very 
small fraction. Most of the particles are so-called cosmic rays consisting of protons, 
helium nuclei and the nuclei of the heavier elements such as carbon, oxygen and iron. 
Just as gamma rays, cosmic rays initiate cascades in the atmosphere. Heavier cosmic 
rays may interact with air nuclei and produce high energy nucleons. High energy pro- 
tons interacting with air nuclei may produce neutral and charged pions. Neutral pions 
have a rather short lifetime and decay, dominantly into photons which may, in turn, pro- 
duce electromagnetic cascades. Charged pions have longer lifetime and may decay into 
muons and neutrinos or may interact with the air nuclei creating secondary high energy 
hadrons and replenish the cascade. Muons, produced in the cascade, may also survive 
to the ground level. The shower stops its growth when secondary high energy hadrons 
and photons can not be produced. 

Because cosmic rays are charged particles they interact with the interstellar and inter- 
planetary magnetic fields and do not provide directional information about their sources. 
Thus, cosmic rays may constitute an unwanted background for a gamma-ray telescope. 
Presence of muons in hadronic cascades is often exploited to differentiate cosmic -ray 
cascades from the gamma-ray ones. 
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2.3 Air shower detection methods. 



The detectors used in high-energy astrophysical experiments are based on those devel- 
oped for laboratory ones. Since the showers extend over large areas large detectors are 
necessary to sample the shower. Cloud/bubble chambers are not suitable for electronic 
data recording and gas-filled discharge tubes are not practical for large area detectors. 
Because charged particles constitute a large fraction of the shower particles scintillation 
and Cherenkov radiation detection techniques are employed in modern air shower de- 
tectors. Cherenkov detectors detect radiation produced when a charged particle moves 
through a dielectric medium at velocity greater than that of light in the medium. Scintil- 
lation counters detect light (luminescence) produced as a result of recombination of the 
electron-hole pairs created by ionizing particles traversing the scintillation medium. 

If the Earth's atmosphere is used as the detection medium this results in the air- 
Cherenkov and "fluorescence" detectors. Air-Cherenkov detectors typically have energy 
thresholds of about several hundreds of GeV, while fluorescence ones can detect high 
energy cosmic rays with energies above 100 (PeV). Such detectors typically have good 
angular resolution but are very narrow field-of-view devices and can operate only on 
cloudless, moonless nights. 

Scintillation arrays have also been built and, due to their sparseness, have rather high 
energy thresholds (typically above several tens of TeV). Such detectors have worse 
angular resolution but can observe the entire overhead sky 24 hours a day regardless of 
weather conditions. 

The goal of the Milagro project is to built a detector sensitive to cosmic gamma rays 
at energies around 1 TeV and capable of continuously monitoring the overhead sky with 
angular resolution of less than 1 degree. 



14 



Chapter 3 

The Milagro Detector 



In a sense it told him nothing that was new, 
but that was part of the attraction. It said 
what he would have said, if it had been pos- 
sible for him to set his scattered thoughts in 
order. 

George Orwell "1984" 

Milagro employs the water Cherenkov detection technique which is widely used in 
particle physics experiments but is new to air shower detection. The use of water as a 
detection medium has several advantages: it is possible to construct a large instrument 
that can detect nearly every relativistic charged shower particle falling within its area by 
observing the Cherenkov radiation the particle produced. At a typical detector altitude, 
there are 4-5 times more photons in an extensive air shower than charged particles. In 
a conventional EAS array these photons are undetected. When these photons enter the 
water, they convert to electron-positron pairs or Compton scatter electrons which, in 
turn, produce Cherenkov radiation that can be detected. Consequently, Milagro has a 
very low energy threshold for an EAS array. 1 

This chapter presents the Milagro detector with its physical and electronic compo- 
nents, event reconstruction methodology and performance characteristics. For a more 
detailed description see references [4] and [3]. 



'Tibet is a conventional EAS array with energy threshold of several TeV. Such a low threshold could 
be achieved only due to its high altitude of 4300 m above sea level [28]. 
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Figure 3.1: Schematic view of the Milagro pond. 

3.1 General description. 

The Milagro experiment is a part of what is now known as the Fenton Hill Observatory 
located at 35.88° North latitude and 106.68° West longitude in the Jemez Mountains 
near Los Alamos, New Mexico. At an altitude of 2650 (m) above the sea level, its at- 
mospheric overburden is about 750 (g/cm 2 ). The Milagro detector, commissioned in 
June of 1999, records about 1700 extensive air shower events per second and is sensitive 
to gamma-showers with energies above 100 (GeV). The duty-cycle of the detector is 
about 90%. The remaining time the detector is off for scheduled maintenance and/or 
when the environmental conditions do not warrant its operation (forest fires, loss of 
electrical power, etc). Milagro is built in a pre-existing 21 (metric) kilo-ton trapezoidal 
water reservoir (see figure 3.1) filled with pure water and instrumented with two hori- 
zontal planar layers of photomultiplier tubes (PMT). The top layer (AS) has 450 PMTs 
arranged on a 2.8 x 2.8 (m) grid, 1.5 (m) below the water surface. The second layer 
(MU) has 273 PMTs located under about 6 (m) of water on an interlaced 2.8 x 2.8 (m) 
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grid. The photo-tube assembly is buoyant with the weight distribution allowing the 
photo-cathode to face upward when the assembly is submerged and anchored to the bot- 
tom of the pond with a Kevlar string. A reflecting conical baffle is installed in each PMT 
assembly to increase the light collection area and block horizontal and upward traveling 
light. The signals from the PMTs are delivered to the data-acquisition (DAQ) system for 
processing and recording. A high-density polypropylene liner and cover are installed to 
ensure water-tight bottom and walls of the pond and light impermeability of the whole 
detector. 

3.1.1 Photomultiplier tube. 

As was mentioned above, the Cherenkov radiation produced in the detector volume is 
detected by photo-multiplier tubes. Unlike conventional electro- vacuum tubes where 
electrons are injected into the tube due to thermal emission from its cathode, the injec- 
tion of electrons (photoelectrons or PE for short) into the photo-tube is caused by light 
via the photo-electronic effect. Due to an externally applied voltage, the electrons travel 
towards the anode of the tube. However, on the way they encounter a dynode chain 
which plays the role of an amplifier. When electrons hit a dynode, secondary electrons 
are emitted which bombard the next dynode on their way to anode. In such a setup, 
enormous amplification can be reached with a relatively short dynode chain. 

The amplification is not the only important parameter of a PMT, the others include: 

Spectral Sensitivity: PMT should be sensitive to the wavelengths produced in the 
Cherenkov radiation. 

Quantum efficiency: The ratio of the number of photoelectrons produced to the num- 
ber of incident on photocathode photons is called quantum efficiency. Ideally it is 
equal to unity. 

Time resolution: Time resolution of a PMT is thought to be limited by fluctuations in 
the photoelectron cascade development, especially on its early stages, especially 
between the photocathode and the first dynode. Lower light intensities generally 
lead to poorer time resolution. 

Pre-pulsing: Pre-pulses on the PMT output are thought to occur when photoelectrons 
are produced by the first dynode, exposed to the incident light. Higher light inten- 
sities generally lead to higher pre-pulsing probability. 
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Late pulsing: Late pulses on the PMT output are thought to occur when all photoelec- 
trons are reflected off the first dynode and re-enter the dynode chain producing a 
PMT pulse later than should have. Higher light intensities generally lead to lower 
late pulsing probability. 

After pulsing: After pulses on the PMT output are thought to occur when residual gas 
molecules in the PMT are being ionized by the photoelectrons. The ions hitting 
the photocathode may cause secondary electron emission which would produce 
a secondary pulse. Higher light intensities generally lead to higher after pulsing 
probability. 

Saturation: Saturation is the effect of decreased PMT amplification for higher intensity 
input. This is caused by inability of dynode chain to accelerate increased numbers 
of secondary electrons to sufficiently high energy. 

After testing several PMT models, the Hamamatsu 8-inch 10-stage R5912SEL was 
selected for this application. It has relatively high quantum efficiency (0.2 — 0.25) at 
wavelengths of 325 — 450 (urn), good timing resolution (2.7 (ns) at 1PE), relatively low 
late/pre/after pulsing rates (about 5%) and relatively long linear response (up to about 
75 PE). For a more detailed description of the PMT characteristics see references [4] 
and [33]. 

3.1.2 PMT pulse model, time over threshold. 

Each PMT should provide information about the intensity of light incident on the PMT 
photocathode and the time when the light was registered. Since the total charge in a 
PMT pulse (number of photoelectrons) is proportional 2 to the incident light intensity, if 
the PMT pulse quickly charges a capacitor which is then slowly discharged via a load 
resistor, the total charge in the PMT pulse can be measured by the capacitor discharge 
time. This suggests that the time spent by a pulse over a preset threshold level is asso- 
ciated with the input light level 3 and is the main assumption of the time-over- threshold 
(ToT) method employed by Milagro. The PMT signal can be digitized with logical 
"one" when the PMT pulse exceeds the discriminator threshold and logical "zero" oth- 
erwise. A time-to-digital converter (TDC) attached to such a digital output will record 
the ToT. The beginning of the logical "one" provides the PMT pulse arrival time (T start ). 

2 Provided that the PMT saturation limit is not reached 

3 In fact, in this model ToT is proportional to the logarithm of the number of PEs. 
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Figure 3.2: Illustration of HiToT, LoToT and edge-train 



This method of measuring PMT pulse charge has several advantages over a conventional 
method when the PMT pulse is sent to an analog-to-digital converter (ADC). ADCs usu- 
ally have narrow dynamic range, and are relatively slow and expensive devices. 

Presence of pre- and after- pulses will distort the PMT pulse and it will not conform 
to the ToT model described above. Since strong pulses are more likely to be distorted, 
two thresholds, high (at the level of about 7 PE) and low (at about 1/4 PE) are in- 
troduced in the Milagro electronics. Large pulses will therefore cross both thresholds 
and the time-over-high-threshold (HiToT) is a much better measure of the pulse charge. 
Two close weak pulses will cross only the low threshold leading to excessively long 
time-over-low-threshold (LoToT), but absence of HiToT will flag such signals. To avoid 
use of two TDCs on a single PMT channel a logical exclusive OR operation is executed 
on LoToT and HiToT digital outputs leading to a train of raising and falling edges cor- 
responding to each threshold crossing (see figure 3.2). If a PMT pulse is weak and only 
low threshold is crossed, the edge train contains only two edges (2-edge pulse), if both 
thresholds are crossed — four edges are recorded (4-edge pulse). Each TDC installed 
in Milagro is capable of recording up to 16 discriminator level crossings. 

The train of edges with their TDC counts constitute the raw PMT signal. 
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Figure 3.3: Block diagram of detector electronics 



3.1.3 The Detector trigger. 

All PMT channels in Milagro were manufactured as uniformly as possible, facilitating 
a simple multiplicity triggering logic. Indeed, as an extensive air shower front hits the 
detector a majority of PMT signals arriving at the outputs of the PMT channels will be in 
close coincidence with each other. The coincidence window was chosen to be 300 (ns). 
If more than 60 PMT signals arrived within the window, a trigger was generated to the 
DAQ system. TDC modules are then read out with look back time of 1.5 (fis) and the 
event is saved. It is desirable to trigger the detector at a low multiplicity requirement 
to lower the detector energy threshold. However, lowering it beyond 60 would increase 
the probability of triggering on muon events which is not the goal of the project. The 
generated trigger was sent to a Global Positioning System clock for absolute event time 
readout. 

The TDC readouts from all PMTs channels and the trigger time constitute the raw 
event data and are sent for further software conditioning and processing. 



3.2 Event Reconstruction. 



The ultimate goal of any high energy gamma-ray telescope project is to study the proper- 
ties of the objects which emitted the particles. This means that the characteristic param- 
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eters of the particles must be defined. Such parameters are: arrival time and direction, 
energy of the particle and its type 4 . (The shower impact parameter on the detector (core 
distance) is also an important parameter, which, however, is not related to the source.) 
As was mentioned in chapter 2, the particles of interest (primary) do not reach the de- 
tector level and disintegrate in the Earth's atmosphere creating extensive air showers of 
secondary particles. These secondary particles can be detected and constitute the ob- 
served event. The process of inferring the characteristics of the primary particle given 
the observed event is called event reconstruction. This is a multi-step process which 
requires deep understanding of the structure of the extensive air showers, detector hard- 
ware limitations and statistical nature of detection itself. 

Currently, the signals from the PMTs in the top layer are used for shower direction 
determination and from the bottom — to distinguish photon and hadron induced air 
showers. 

3.2.1 Pre-processing. 

As was described above, the raw event data contains "edge-train information" registered 
by each PMT's TDC in the event. These data are not immediately suitable for primary 
particle characterization because the data is tainted by noise and systematic effects in 
the detector. Systematic effects include systematic off-sets of TDCs on different PMT 
channels (called time pedestals), TDC conversion factors (number of TDC counts per 
unit time) and delayed electronics response to lower PMT signals compared to higher 
ones (called electronic slewing). These systematic effects are studied with the help of 
the calibration system (see appendix B) and can be taken into account. Noise effects are 
random by nature and thus are more difficult to study and correct. Noise could be due to 
signals not associated with the main shower event (thermal/electronic/radioactive noise, 
non-shower particle hitting a PMT) or partially recorded edge-trains due to hardware 
constraints. 

Noise filtering. 

The main purpose of the edge-finding filter is to check that each PMT signal in an event 
conforms to the PMT signal model described in section 3.1.2. This should eliminate 
some thermal/electronic noise and partially recorded signals. The behavior of all PMT- 

4 major types are photon and hadron 
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electronic channels was studied in great detail and based on that, a set of criteria was 
developed which would select viable signals. (See [27].) 

If a PMT signal does not satisfy the criteria, an attempt is made to convert the signal 
to the proper form. This is done by checking the number of edges, their polarity 5 and 
timing within the PMT signal. This filter is applied for each PMT in each event regis- 
tered by the detector. After such filtering only about eight percent of all PMT signals 
are considered as unrecoverable and are discarded. 

A completely different problem arises when valid PMT signals from non-shower 
particles are recorded with the main shower event. Presence of these signals will degrade 
the quality of event reconstruction as such signals do not carry any useful information 
about the shower. An idea of a method for filtration of such PMT signals was first 
proposed in [26] and [50] and then used in [23] and is based on the fact that PMT 
signals produced by a shower must be causally related, i.e. the time interval between any 
two PMT pulses multiplied by the speed of light in water should not be larger than the 
spatial distance between the PMTs. If a PMT signal is causally disconnected from the 
main shower event, it should be discarded. Unfortunately, at the moment of writing, this 
idea is not developed enough to be a part of the standard Milagro event reconstruction. 
This filter is applied to calibrated event. 

TDC conversion factors. 

TDC conversion factors were monitored with the help of the calibration system (see 
appendix B) and were found to be stable. Of most importance is the fact that all TDC 
modules operated at a common conversion rate of 2 counts per nanosecond with ex- 
tremely high precision (see section B.2.1). This means that TDC counts can be used as 
time measure directly and there is no need to convert TDC counts to time for each PMT 
channel separately. This simplified the structure of the reconstruction code. 

Raw to Calibrated event. 

As was mentioned before, time response of PMT-electronic channels is dependent on 
the light intensity input. Since each PMT signal remaining after filtration conforms to 
the signal model, it is possible to correct for the effects of electronic delays to signals 
of varying strength. Based on that and auxiliary data obtained from the calibration 
process (see appendix B) the measured Time-over- Threshold (LoToT for weak signals 

5 polarity means correct sequence of rising/falling edges 
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and HiToT for strong ones) can be converted to PMT pulse arrival time (T start ) and 
number of photoelectrons emitted from PMT photocathode (PE). PMT coordinates and 
observed T start and PE for each PMT in an event constitute a calibrated event and contain 
all information needed for event reconstruction. 

3.2.2 Processing: angle, time, energy, type. 

Even though at this processing stage all information obtainable from the PMTs is known, 
to reconstruct the particle characteristics the general structure of the showers and detec- 
tor capabilities should be taken into account. For instance, since the PMT efficiency is 
only about 20 percent, there is no guarantee that the observed PMT signal is generated 
by the shower-front particles. Particles trailing the front may generate a PMT pulse too, 
but if the PMT happened to register the shower-front particles, the PMT pulse might 
be discarded as the ToT pulse model does not allow for more than a single PMT pulse 
in a shower event. Another example is that it is almost impossible to differentiate a 
low energy shower with small detector impact parameter from a high energy one with 
large impact parameter without knowledge of the shower structure. Thus, any method 
of event reconstruction must take into account detector and shower features. 

Shower impact parameter. 

As discussed in chapter 2 while the primary particle impact parameter does not provide 
any information about the source and the particle it created, it helps to understand the 
detector response to the shower produced. 

Currently, the PMT PE distribution in the top layer in an event is analyzed to infer 
the location of the shower core. If the decision is made that the core is inside the Milagro 
pond, a PE-weighted average of PMT positions is used as the core location, while if it 
is decided that the core is outside the pond, it is placed at the distance of 50 (m) 6 from 
the center of the pond. The direction to the core, in the later case, is reconstructed by 
connecting the center of the pond with the V P_E-weighted PMT positions. The decision 
of whether the core is inside or outside the pound is made based on the radial profile 
distribution of the number of PEs observed in the top layer PMTs. 

The information inferred about the shower core is used in the sampling correction, 
angular and energy reconstruction. Full details of this method are described in [48]. 

6 Computer simulations indicate that this is the most probable core distance for the showers which 
trigger the detector and have cores outside of the detector. 
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Sampling correction. 

A great care has been take to eliminate systematic and random effect in the detector on 
event reconstruction. There is, however, a remaining one. This has to do with the finite 
probability of a PMT-electronic channel to detect light. Thus, the light, produced by 
the shower particles may be lost. The situation is complicated by the fact that showers 
have thickness and detection of trailing particles, if interpreted as the shower front, will 
degrade the quality of the angular reconstruction. 

Luckily, knowing that the thickness of the shower is a function of impact parameter 
and that the number of particles in the shower falls off quickly with longitudinal distance 
from the shower front, the amount of light produced by the trailing particles is generally 
lower than by the front of the shower. Using that knowledge, the shower sampling effect 
can be observed based on measured light level at a given PMT and its distance from 
the shower axis and PMT pulse time can then be corrected to represent the shower front 
arrival time. 

The Milagro sampling correction has been developed based only on number of PEs 
registered by a PMT in [35] and [22] and assumes that the shower arrived vertically 
when the impact parameter is equal to the core distance. 

Time of event. 

The time of an event is recorded as time of arrival of the PMT multiplicity trigger and is 
read from a GPS clock. 

Angular reconstruction. 

After detector sampling effects have been taken into account, the obtained PMTs' T start 
times represent the best knowledge of the shower front. Knowing that the shower front 
forms a paraboloid, its main axis can be found and will give the arrival direction of the 
progenitor particle. 

The algorithm utilized by Milagro first assumes that the shower arrived vertically and 
given the shower core position the curvature of the shower front can be "taken out" with 
what is called "curvature correction". 7 Following that, the shower direction is sought as 
the directrix of the plane fitted to the PMTs' T start times ("time-lag" method) using a 
weighted x 2 -method. (See for instance [15].) The weights for the x 2 -fit are prescribed 

7 The functional form of the curvature correction was obtained from data and computer simulations 
in [35]. 
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based on the number of PEs observed, as the quality of PMT time resolution increases 
with increase in the input light level [35]. 

Energy reconstruction. 

Energy estimation is based on the amount of light deposited in the detector, distance to 
the core and the angle of the shower arrival and relies heavily on computer simulations 
of the shower propagation in the atmosphere and in the detector. At this time, primary 
particle energy is not being inferred in online data processing. 

Primary particle type identification. 

Because of their hadronic cores, air-showers generated by incident cosmic rays develop 
differently from purely electromagnetic cascades. The probability of photons to produce 
electron-positron pairs is several orders of magnitude higher than that of any process that 
might lead to muon production. In contrast, interactions of high energy hadrons with 
atmospheric nuclei lead to the production of charged pions which may decay into muons. 
In addition, multi-GeV hadronic particles may also survive to the ground. Simulations 
indicate that 80% of proton and only 6% of photon induced air showers that trigger 
Milagro will have at least one muon or hadron entering the pond. 

Hadrons that reach the ground level and produce hadronic cascades in the detector 
or muons that penetrate to the bottom layer will illuminate a relatively small number 
of neighboring PMTs in that layer. Photon induced showers, on the other hand, gen- 
erally will produce rather smooth light intensity distributions. Based on this simple 
observation a technique for identification of photon/hadron initiated showers has been 
formulated [51, 5] and according to computer simulations can correctly select about 
90% of hadron initiated showers and about 50% of photon induced ones. 

In a search for sources of high energy photons where hadron initiated showers rep- 
resent unwanted background, the proposed identification scheme will allow increase of 
signal to noise ratio. 

3.2.3 Post-Processing: Analysis Techniques. 

After characteristics of the primary particle have been established further analysis has 
to be done, based on the concrete task under investigation. While many different tasks 
use similar techniques to answer stated questions, many employ unique methods. For 
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this reason the discussion of methods and algorithms used in the present work is delayed 
until chapter 4. 

3.3 Detector performance and simulations. 

After a device have been built and tuned, it is desirable to test its operation and gauge 
its response. Usually, this is done by comparing the device's response with the expected 
one given a known input signal. Needless to say that such a test is not possible to 
perform with Milagro due to unavailability of controllable test sources of high energy 
particles above the Earth's atmosphere, and one is forced to resort to computer simula- 
tions to estimate the detector performance. A simulated extensive air shower is sent to 
a simulated detector. The output of the simulated detector is sent for standard analysis 
and the result is compared with the input primary particle parameters. 

The air shower simulation is done with the CORSIKA package [30] in the standard 
US atmosphere down to the detector level. The simulated shower front is then input 
into the GEANT-based detector simulation package. The output of this procedure is the 
Milagro "calibrated" event which can be sent for the standard particle characteristics 
reconstruction described earlier. 

The most important parameters of the detector which are obtained based on com- 
puter simulations are angular resolution, energy response, impact parameter informa- 
tion, particle type identification quality and the detector's effective area. 

Extensive air showers were generated over an energy range of 100 (GeV) to 
100 (TeV), with zenith angles ranging from to 45 degrees and core locations uni- 
formly distributed over 1000 (to) radius around the detector. Probability of triggering 
on a shower with energy outside the selected range or with higher zenith angles is very 
small which motivated the choice. 

Angular resolution is characterized by the difference between the reconstructed and 
the known input particle direction. The overall accuracy of the angular reconstruction 
is believed to be 0.75°. The report [52] suggests that the energy of an incident particle 
can be reconstructed by Milagro with a fractional error of about 50% for particles with 
energies above 1 (TeV). 8 Core location is reconstructed with error of about 20 (to) if 
the shower lands on the detector and with error of about 50 (to) otherwise. 

8 The same report also implies that since the quality of energy reconstruction relies on the quality of 
the core reconstruction, it is not possible to reach the 50% energy resolution with the current Milagro 
configuration. An upgrade with an outrigger array is necessary. 
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3.3.1 Effective area. 



As was already mentioned, a shower event can be detected even if its core lands outside 
the detector. This leads to the notion of effective area as the area of imaginary detector 
which has perfect sensitivity to events which land on it and zero outside. This parameter 
describes sensitivity of the detector to particles of different type, energy and arrival 
direction. 

If N (k, E, 9) showers induced by particles of type k are simulated with core lo- 
cations uniformly distributed over sufficiently large area A , local arrival directions 
(6, 6 + SO) and energies in the range of (E, E + 5E) then the effective area A k (E, 0) 
can be computed using the number of events N t (k, E, 6) which satisfy detector trigger 
condition in the simulations: 

A (F Cj\ N t (k,E,e) 

ME > 0) - N (k,E,e) Ao 

Base simulations of proton and photon initiated showers A k (E,9) were obtained 
where 9 is the local zenith angle only. 
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Chapter 4 

Analysis Techniques 



...if all others accepted the lie which the 
Party imposed — if all records told the same 
tale — then the lie passed into history and be- 
came truth. 

George Orwell "1984" 



4.1 Coordinates on the Celestial Sphere. 
4.1.1 Celestial Sphere. 

Because the stars are distant objects, they appear to lie on a sphere concentric with 
the Earth. This imaginary sphere is known as the Celestial Sphere. Astronomy uses a 
number of different coordinate systems to specify the positions of celestial objects and 
only those relevant to this work ones are discussed here. 

The Celestial sphere has North and South Celestial poles as well as the celestial 
equator which are projected reference points of the same positions on the Earth's sur- 
face. A coordinate system which is based on these reference points on the celestial 
sphere is called the equatorial celestial coordinate system and is similar to the geo- 
graphical coordinate system on the Earth's surface. A point on the celestial sphere can 
be described by two coordinates named "declination" (5) and "right ascension" (a). The 
declination of a star is the analog of the latitude and is the angular distance from the star 
to the celestial equator. Right ascension is the analog of longitude with the zero of 
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Figure 4.1: Definitions of main points and arcs on the Celestial sphere 



right ascension at the point of vernal equinox. 1 Because of the Earth's rotation, right 
ascension and declination are not measurable directly by the ground-based observer and 
additional local to the observer coordinate system is introduced. The local coordinates 
are azimuth (A) and zenith distance (z) which can be converted to declination (5) and 
hour angle (H). The list below defines the main points and arcs on the celestial sphere 
which are illustrated on figure 4.1. 

C — Observer 

CP — Axis parallel to the axis of rotation of the Earth passing through the observer C. 
P, P' — North and South Celestial poles. 

Z — Zenith, CZ is the continuation of the plumb line at observer C. 

Horizon — intersection of a plane perpendicular to CZ at point C and the celestial 
sphere. 

'Vernal equinox is the point where the Sun crosses the celestial equator on its south to north path 
through the sky and is stationary in space. 
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Local Reference Meridian — The arc PZP' of the great circle 2 containing points P, 
Z and P'. 

N,S — North and South on the horizon as defined by the intersection of the great circle 
PZP' with the horizon. 

PZ — LPCZ = 7r/2 — 4>, cj) is geographical latitude of the observer on the Earth. 

Celestial Equator — intersection of the plane perpendicular to CP at point C with the 
Celestial sphere. 

X — A celestial object on the sky (a star). 

ZX — Zenith distance (z) of the star X is the angle LZCX. 

LPZX — Azimuth (A) of the star X is the dihedral 3 angle between the reference merid- 
ian and the ZCX plane measured from North towards East. 

PX — IPCX = 7r/2 — 5, declination (5) of the star is the angle between CX and the 
Celestial equator. 

LZPX — Hour angle (H) of the star is the dihedral angle between the reference merid- 
ian PZ and the PCX plane measured from South towards West. 

T — Point of vernal equinox 

PTP' — Celestial reference meridian. 

TPX — Right ascension (a) of the star is the dihedral angle between the TCP and 
XCP planes. 

Given the definitions above, the law of cosines for the trihedron 4 applied to the 
spherical triangle ZPX two times yields the relationship between the (A, z) and {8, H) 
coordinate systems: 

2 A great circle is a section of a sphere that contains a diameter of the sphere. 

3 The dihedral angle is the angle between two planes and is defined as the angle between their normal 
vectors. 

4 Three vectors with common vertex, often called a trihedral angle since they define three planes. 
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cos(7r/2 — 5) = cos z cos(7r/2 — 0) + sinzsin(7r/2 — 0) cos(27r — A) 
cos £ = cos(7r/2 — 5) cos(7r/2 — 0) + sin(7r/2 — 5) sin(7r/2 — 0) cosH 

{sin 5 = sin cos z + cos sin z cos A 
tan// = - sin z sin A 

cos <p cos z—smz sin cos A 

Since the local reference meridian is defined relative to the Earth, due to Earth's 
rotation the hour angle of a fixed in space point will grow with time (that is why it is 
called hour angle) while the local coordinate declination will remain constant. The hour 
angle of vernal equinox H r links the local observer's coordinate system (H,5) and the 
celestial equatorial coordinate system (a,S) by providing the position of a fixed point 
(vernal equinox) on the celestial sphere in local coordinates: a = H T — H. Hour angle 
of vernal equinox is also called the local sidereal time since it should be consistent with 
the observer's geographical longitude and the time required for one Earth's revolution, 
called a sidereal day. In contrast, the solar (or universal) day is defined as time between 
two consecutive appearances of the Sun on the local reference meridian. The solar day 
is longer than the sidereal one due to Earth's rotation and orbital motion around the Sun, 
though both days are divided into 24 hours. 

4.1.2 J2000 reference. 

Because the Earth's rotation is not uniform, its axis of rotation is not fixed in space and 
even its shape and relative positions on its surface are not fixed; because the introduced 
celestial equatorial coordinate system follows the motion of the Earth's pole and equa- 
tor, the coordinate grid "drifts" on the surface of the celestial sphere. 5 Therefore, the 
introduced coordinate system provides only apparent right ascension and declination of 
the stars at the observation moment. 

To solve this problem, all coordinates on the celestial sphere are reported relative 
to the position of the Earth's pole and equator at specified moments of time which are 
called epochs. Each epoch lasts for 50 years and the current one is defined with respect 
to the Earth's position at noon on the January 1, 2000. Thus, the apparent celestial 
coordinates need to be reduced to the J2000 reference. 6 

5 These drifts include, but not limited to precession, nutation, celestial pole offset and polar motion. 
6 The major contribution to the "drift" of a celestial reference frame is due to the Earth's pole preces- 
sion. Newcomb (Newcomb, S. Astron.J. 17, 20 1897) derived the formulae for precession parameters 
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4.1.3 Diurnal parallax. 



The Equatorial coordinate system had been defined under the assumption that the ob- 
server is located at its origin — the center of the Earth. All observing stations, however, 
are located on the Earth's surface. Due to Earth's rotation, the observing station moves 
and the observation of a celestial body is being made from different points in space. 
This will cause an apparent difference in position of celestial body when made at dif- 
ferent moments of time. The effect is called diurnal parallax. For measurements of 
distant stars this has a negligible effect, but there could be a substantial diurnal parallax 
on objects inside the Solar system. Diurnal parallax on the Moon, for example, can be 
as large as 0.95°. 

4.1.4 Milagro event coordinates. 

The local hour angle and declination of an event on the celestial sphere are calculated 
from the zenith and azimuth which are provided by the event reconstruction section 3.2.2 
(see also [19] for a discussion on local coordinates). Local sidereal time as well as 
the geographic coordinates of the detector can be obtained from a Global Positioning 

Ca(*)> ZA(t) and which specify the position of mean equinox and the equator of a date with respect 
to the mean equinox and equator of the initial epoch. Astronomical Almanac on page B18 supplies these 
parameters for the J2000.0 epoch in degrees: 



where T stands for the time from the basic epoch J2000.0 in Julian centuries, T — (Julian Day — 

2451545.0)/36525. 

If subscript refers to the coordinates at the epoch J2000.0 and no subscript to the epoch of the date, 
the transformation formulae are: 



( A = 0.6406161T+ 0.0000839T 2 + 0.0000050T 3 
z A = 0.6406161T+ 0.0003041T 2 + 0.0000051T 3 
9 A = 0.5567530T- 0.0001185T 2 - 0.0000116T 3 



sin(a — za) cos S 
cos(a — za) cos S 
sin S 



sin(a + CO cos<5 

cos(ao + CO cos 9a cos So — sin 9a sin So 
cos(ao + CO sin 9a cos So + cos 9a sin <5q 



sin(a + (a) cosS 
cos ( a + CO cos<5 
sinJo 



sin(a — za) cos S 

cos (a — za) cos 9a cos S + sin 9a sin 5 
— cos(a — za) sin 9a cos 5 + cos 9a sin S 
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System receiver which facilitates the conversion from local to celestial coordinates. In 
Milagro, the coordinates of reconstructed events are reduced from the epoch of date to 
the J2000 reference in real time and are saved to disk for further processing. 

As will be clarified in the sections to follow, the signal processing method employed 
in this work expects the event coordinates in a local reference frame. Thus, even though 
the events are saved in J2000 reference which seems to be convenient, the conversion to 
J2000 must be undone during the offline/online signal processing. 7 



In counting-type astrophysical experiments, the brightness of a particular point in the 
sky is characterized by the number of events observed from that point during the ex- 
posure time. Such experiments measure the density of events on the surface of the 
Celestial sphere. Therefore, the procedure used to generate the sky images (projections 
of a sphere onto a plane surface) should conserve the density of events. To meet this 
requirement an equal area projection of a sphere onto a plane must be used. This re- 
quirement, however, is not enough to uniquely fix the projection, and several different 
projections are available. It is crucial to understand that any area preserving projection 
is not conformal and might distort the distances and/or directions on the map. That is 
why different area preserving mappings should be used for different tasks. For example, 
the sinusoidal projection 8 has minimal distortions near the equator and that is why is it 
very convenient for galactic plane studies. When the same mapping is used for an object 
far from the galactic equator, linear distortions become significant. 

7 It would be prudent to save the local hour angle and declination of the registered events during 
the online realtime processing. This would force the coordinates of celestial bodies which are known in 
J2000 from catalogues, to be reduced to epoch of date then to be reduced to the apparent Right Ascension- 
Declination by application of the parallax correction (if necessary) then to local hour angle-declination 
using local time. This approach would save some computer time during online and offline data processing 
because the Milagro event rate is above 1 (kHz) and the detector angular of resolution (several tenths 
of a degree) allows for rare (once per several seconds (in 24 seconds the Earth rotates on 0.1° of arc)) 
computation of local coordinates of the celestial bodies. 

8 The Sinusoidal Equal Area Projection is defined as: 



4.2 Sky Mapping. 




where (I, b) are galactic coordinates 
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Figure 4.2: Concept of the auxiliary coordinate system on the Celestial sphere cen- 
tered on L (left) and corresponding sky projection (right). The "y"-axis of the sky 
projection always points to a pole M and the circles are the lines of x — const at 



In addition to the previously discussed constraints, all events with identical spatial 
orientation with respect to the point of interest must be mapped into a unique location. 
This is especially important if the point of interest moves on the Celestial sphere. How- 
ever, a simple algebraic difference in Celestial coordinates of any two points does not 
define their relative spatial orientation. One way to address this problem is to introduce 
auxiliary coordinates on the Celestial sphere: an analog of latitude (x) and longitude (£) 
which are measured with respect to the preselected point. (See figure 4.2 and appendix C 
for the definition of the (% — £) coordinate system.) Using these coordinates, the sky 
image centered on the selected point can be produced with the help of the Azimuthal 
Equal Area Projection in the polar case defined as: 



This mapping not only satisfies the above requirements, it has several other impor- 
tant features such as conservation of the directions as seen from its origin, the locus 
of points equidistant from the center of the mapping is projected into a circle (see fig- 
ure 4.2) and it can be easily oriented along the lines of the Earth's magnetic field. 



X = t/6, tt/3, tt/2. 




(4.1) 
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4.3 Statistical Nature of Signal Establishment. 



In a typical counting type astrophysical experiment during observation time ti the num- 
ber of events observed due to some physical process is N\. Assuming that an event 
contains no information about any other one, the number of observed events is a random 
variable which is distributed according to the Poisson distribution. 9 In other words, the 
probability to observe exactly N\ events during time t\ is given by: 

where A has a meaning of average event rate. 

However, an observed event could be due to either a source or background. Since 
the average count rate due to background is not known, based on this one observation 
it is not possible to decide whether there were any "source" events observed. Therefore 
without altering the conditions of the experiment a second measurement N 2 during t 2 
is made where it is believed that all observed events are due to background only. Now, 
a decision should be made as to whether there is a difference between these numbers 
which can be interpreted as a detection of a source. Since the observed JVi and iV 2 are 
random numbers, this question should be approached from the statistical point of view. 
Note that a statistical test cannot verify that a given hypothesis 10 is true or false, but can 
only suggest which of the two or more hypotheses is the more plausible explanation of 
the observation. 



4.3.1 General test construction. 

Suppose 11 that a result of an observation is described by the values of n variables: 

\x) = x±, x 2 , ■ ■ ■ , x n 

The {x} may represent outcomes of n repeated measurements made under identical 
conditions or a sample from a population. All possible outcomes of a measurement are 



9 Some of the properties of the Poisson distribution are discussed in appendix A 
10 Any statement concerning the unknown distribution of a random variable is called a statistical hy- 
pothesis. 

"This subsection is based on the section II of the paper [40]. 
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said to form a sample space W. A hypothesis about the origin of the observed events H 
defines the probability of occurrence of every possible observation 



p(x 1 ,x 2 , ...,x n ) 



and thus, the probability that the observed event will fall into some region w of all 
possible outcomes is 



Of course, P{W) = 1. Different hypotheses H with their corresponding probability 
distributions p({x}) will be endowed by the same subscripts, such as H and p . 

A statistical test is formulated so that all prior knowledge strongly supports H called 
the null hypothesis. Hypothesis H is rejected if the observed event {x} lies within a 
certain critical region w c and accepted or doubted otherwise. Such a formulation of a 
test implies that it is possible to reject H when, in fact, it is true. The danger of falsely 
rejecting the null hypothesis is characterized by the error of the first kind or significance 
£ c and: 



The choice of the value of £ c depends on the penalty for making the error, therefore, 
the risk £ c must be set in advance, not after the results of a measurement are available. 
Even though the error of the first kind can be chosen to be arbitrary small, the equation 
£c = Po(w c ) has, in general, infinitely many solutions on configuration w c with the same 
level of significance £ c . 

Since H is being tested, it implies the existence of an alternative hypothesis Hi or 
there would be no question about H . n But, as the risk £ c is required to be smaller and 
smaller, the risk ( of accepting H , when Hi is true may increase. This error is called 
the error of the second kind and is given by: 



The two errors £ and ( can rarely be eliminated, and in some cases it is more impor- 



While it may not be constructive, "Hq is false" is an admissible alternative hypothesis. 




n 




n 




n 
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tant to avoid the first, in others — the second. When H and Hi are specified it is the 
choice of the critical region which allows control of the errors. 

A prescription to resolve the apparent vagueness in the provided formulation of the 
test was proposed in [40]. It is proposed that given the two hypotheses H and Hi and 
the desired risk level £ c , the corresponding best critical region w b c est minimizes the error 



If based on the outcome of the experiment, the observed x is inside of the critical 
region w b c est , it is said that the null hypothesis is rejected in favor of the alternative one 
with significance £ c and power (1 — £). If, however, x w best , it is said that the null 
hypothesis is not rejected in favor of the alternative one with significance £ c and power 



Often, however, rather than use the full data sample {x} it is convenient to define 
a test statistic 13 U. Each hypothesis for the distribution of {x} will determine a distri- 
bution for U, and a specific range of values of U will be mapped to a critical region in 
W-space. In constructing U one attempts to reduce the volume of data without loss of 
the ability to discriminate between different hypotheses. 

4.3.2 Testing a composite hypothesis. 

If the hypothesis being tested does not specify the probability of occurrence of every 
possible observation, it is called a composite hypothesis. It will be assumed that the 
composite hypothesis depends on an unspecified parameter A as: 



As before, the null hypothesis should be rejected if the observed event lies within a 
critical region w c . 

In order to control the error of the first kind £, the critical region must satisfy: 



for every value of the parameter A . In other words, the error of the first kind should 
not depend on the unknown value of the parameter Ao- If such critical regions exist, it 
is necessary to choose the best one which minimizes the error of the second kind. It 

13 Statistic is a random variable which is a function of the observed sample of data. 



C- 



(i-C). 



p(xi,x 2 , A) 




n 
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should be noted that the error of the second kind may, in general, depend on the values 
A and Ai of the alternative hypothesis pi(x± : . . . , x n ; Ai). 

This problem has been solved in [40] for a special class of the null hypotheses when 
Po({x}; A ) is infinitely differentiable function of A in every point {x} £ W and the 
function po({x}; A ) satisfies the equation: 

# = A + 50, = ^MMiM (4 .2) 
aAo aAo 

and the coefficients A and B are functions of A only and do not depend on {x}. It is 
shown in [40] that the best critical region w b ^ st is constructed of pieces of hypersurfaces 
= C = const such that: 

g #TTT >^ VWg»« (4-3) 
Po({z};A ) 

where q is a constant whose value is governed by w b c est chosen subject to constraint: 

ic / Po({x}; \ )dx = / Po({x}; \ )dx (4.4) 

4.3.3 Significance of a measurement. 

In as much as an attempt is being made to identify the presence of a source, the null- 
hypothesis H will be formulated in the following way: 

The source is not present. The results Ni and N 2 of two independent 
observations come from a single Poisson distribution with parameter A. 

with an alternative hypothesis Hi that: 

The independent counts N\ and N 2 come from Poisson distributions 
with different parameters Ai and A, correspondingly. 

Mathematically, if H is true the probability po(N 1 , N 2 ; A) to observe N± and N 2 is: 
while, if Hi is true the probability pi(Ni, N 2 ; Ai, A) is: 
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p lW , J V 2 ;A 1 ,A)=(Mp e -^(^e-«. 

where the values of A and Ai are unspecified and the only requirement is that Ai 7^ A. 

The formulated H satisfies the conditions of a theorem presented in [40] which 
states that there exists the best critical region w^ est corresponding to significance £ c 
independent of the value of the parameter A. 

Following the algorithm for construction of the best critical region from [40], the 
equations (4.2) and (4.3) become: 



Pi 
Po 



> q 



Ai 



(f) = N t = Nx + N 2 = const 

\ iVi 



e -(Ai-A)ti > q ^ 



Ni > iV 6 Ai > A 
Ni < Nt., Ai < A 



It is thus clear that the best critical region for testing Ai = A against Ai 7^ A does not 
exist, however, it does exist for testing Ai = A against Ai > A or Ai < A separately. 

The value N% corresponding to the error of the first kind £ is found as the solution of 
the equation (4.4): 



£ HkloPo(k, N t — k;X) — Etm Po(k, N t -k;X), X, > A 



eEfi Po(fc, N t — k;X) — E2oPo(k, N t - k; A), X 1 < A 



Immediately, it should be noted that the solution N% does not depend on the values 
of the parameters A x and A and the best critical region exists for the H with regard to 
all alternative hypotheses Hi. After explicitly writing the probability po(N 1: N 2 ), one 
arrives to the following equation on Nf 



i = {I + a)- N ^t N C k Nt a\ Ai > A 



H = (l + a)- N <EkloC k Nt a k , Ai < A 



a 



h/t 2 >o, 



m\(n — m)\ 



The error of the second kind ( can be computed as: 
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C = ££ t=0 Efio "Piik, N t - k; Ax, A), X 1 > A 
C = E^ =0 Eki Nl+ iPi(k,N t - fc; Ai, A), Ai < A 



C = e 



— Aiti — A<2 



iv c -i (^) fc 



Ai > A 



'fc=0 fc!(JV t -fc)P 



— Aiti — A<2 



Chilli*: 
1 At 2 -> 



Ai < A 



■N s +1 fc!(jV t -fc)P 



The explicit solution for the critical region is needed if an ability to compute the 
error of the second kind ( is desired. As expected, this error will depend on the values 
of the parameters A and A x of the alternative hypothesis Hi. It is, however, possible 
to decide if the null hypothesis should be rejected or not without the explicit solution. 
To do this, iVg must be set to Ni and £ must be computed from the equations above 
using N t = (Ni + N 2 ). If it is found that the £ obtained in this fashion is smaller than 
the critical value £ c , the null hypothesis should be rejected and should not be rejected 
otherwise. 

As will be explained below, because the procedure for setting an upper limit is based 
on the error of the second kind (, expression for which is not known in a closed form, a 
"practical" statistic which was proposed in [34] is considered in this work: 



The denominator in (4.5) is the maximum likelihood estimate on dispersion of (Ni — 
aN 2 ) given the null hypothesis is true. Then, under the null hypothesis the mean value 
of the statistic U is zero and the dispersion is equal to unity. If both Ny and aN 2 have 
not deviated far from the expected value of Xti, then N x and aN 2 can be regarded as 
coming from Gaussian distributions with the means equal to Xti and dispersions Ati 
and aXti correspondingly (See discussion on Gaussian limit to Poisson distribution in 
appendix A.). Hence, the values of the statistic U are distributed according to a Gaussian 

u 2 

distribution with zero mean and unit variance po(u) = -k=e ~. This statement is valid 



U 



Ni - aN 2 



a = ti/t 2 > 



(4.5) 



fi^Ni + N 2 ) 
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for all u < Uq U : 



\u\ < \uq\ « min ^3Qa(l + a) 2 (N 1 + N 2 ), ^36«- 3 (l + a) 2 (A^ + JV 2 )) (4.6) 

If, however, the i?i(Ai) is true, the U will have approximately Gaussian distribution 
with unit dispersion and shifted mean: 



where u\{\\) is monotonically increasing function of Ai and is equal to the average 
value of U computed when Hi(Xi) is true. 

/\ \ (Ai — Ajtx Ai — A r— 



Let us define the critical range of values of the statistic U corresponding to signifi- 
cance £ c in the following way (see figure 4.3 for an illustration): 

If Ai > A: u> u c , £ c = /+°° p (u)du. 

If Ai < A: u< u c , £ c = J^poi^du. 

For the reasons of tradition, in astrophysics, it is customary to report the level of sig- 
nificance not as probability £ c , but as "number of sigmas" u c which motivated the choice 
of statistic. The table 4.1 provides the translation between the critical value u c and the 

2 

significance £ c with the approximation error on £ not exceeding ^== J^ 00 e^^du. 



4.3.4 Setting an upper limit. 

Some times, when the null hypothesis can not be rejected based on the results of a test 
and there are several alternative hypotheses available, it is meaningful to ask the ques- 
tion of which of the alternatives provides error of the second kind larger than ( u . For in- 
stance, in the case considered here, there are many alternative hypothesis parametrized 

14 The value of uq is obtained by substituting k with N1.2 and A with corresponding maximum likeli- 
hood estimates ^rjr^ i 2 into the equation (A.2). 
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Table 4.1: Significance £ c and corresponding critical value u c . 




u c «i(Ai) 



Figure 4.3: Critical region illustration for the statistic U when Ai > A. 
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by the values of A! 15 . Since each alternative hypothesis if^Ai) defines a probability 
distribution p±(u, Ai) in the sample space, the error C(^i) is ( see figure 4.3 for an illus- 
tration): 



For the case of Ai > A, C(^i) is monotonically decreasing function of Ai. Therefore, 
A" corresponding to the largest allowed error ( u is called the upper limit on Ai (C(A") = 
( u ). It means that the probability of making the error of the second kind by accepting 
the null hypothesis when in fact one of the alternative hypotheses with parameter Ai 
(Ai > A" > A) is true is less than ( u . For a discussion on the upper limit construction 
procedure, please, see appendix E. 

4.4 Background estimation. 

One method of searching for a source is by counting the number iV on (0) of events 
from local direction (0; + d<d) while it was exposed to a source region f2 on the sky 
(the "on-source" bin) and comparing it with the number of background events A^ n (0) 
expected from this region. The number of background events expected from the "on- 
source" region is given by: 



where -R&(0, t) is background event rate from local direction at an observation mo- 
ment t and 



Since the function Rb(Q, t) is not known a priori, it should be determined from the 
observed data. To accomplish that one is forced to introduce some assumptions about 
the structure of /?{,(©, t). Probably, the most natural simplification comes from the as- 
sumptions that the background events are distributed uniformly on the sky (their distri- 

15 Remember that Ai is not a source strength, but merely the average event count rate due to possible 
presence of a source. 



IfAi> A: C(Ai) = J^oPi(«,Ai)d«. 



If Ai < A: C(Ai) = J^°Pi(«,A 1 )d«. 
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bution is independent of local coordinates 0) 16 and that the conditions of the experiment 
(hardware, software, field of view, everything) remain constant (at least for periods of 
time long enough to allow measurement of R b (Q, t) to the necessary accuracy). Then, 
Rb(Q, t) can be factorized: 

Rb(&,t) = G(Q)-R b (t) 

where Rb(t) is overall event rate due to background only and is independent of local 
coordinates and G(Q) is the detection efficiency of the instrument and does not depend 
on time. Thus, the problem of determining R b (Q, t) is reduced to the one of R b (t) and 
G(0). 

Knowing that the number of background events expected from any point on the sky 
at some time is: 

dN b (Q,t) = R b (Q,t)d&dt = R b (t)G(&)d&dt 

the total number of background events A^ ut (©) which are to be observed within some 
large time T from the local direction outside of the source region 17 18 and the total 
rate -R^(t) from all viewed sky except for the source region are: 

f ;vl(0) = f T <j>(e,t)G(e)R b (t)dt = G(e)j T d>(e,t)R b (t)dt 

\ R b out (t) = 1 0(0, t)G(Q)R b (t)de = R b (t) j 0(0, t)G(e)dQ 

The functions N b ut (Q) and R b out (t) are not distorted by the presence of any source 
(by assumption) and can be measured directly. Then, the set of equations can be solved 
for R b (t) and G(Q) numerically with the initial approximation to R b (t) taken from the 
observed total event rate. 

Thus, the expected number of background events in the source region can be found 
from: 

16 Charged particles which form the background are isotropized by galactic and inter-galactic magnetic 
fields. 

17 The outside region should not contain any known source in the field of view of the detector. The 
events from other sources and their source regions should be removed from the analysis as they would 
bias the background estimation. 

18 Due to the Earth's rotation, the local direction 9 will fall within the source region £1 at some times 
and outside at the others. If for a particular source region 0, 9 is exposed to only, N^ n can not be 
estimated with the presented method. 
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N b on (Q)= f [1 -<!>(&', t)}R b (t)G(Q')dQ'dt 

4.5 Performing test for a source presence. 

If N on (Q) and N b n (Q) are the number of events observed from the local direction 9 
inside of some bin f2 in the on-source and off-source measurements respectively, the 
value of the statistic U from (4.5) is: 



v /ivi(e) + a(e)iv on (e) 

where 

N 1 = N on (Q) aN 2 = N b on (e) «(©) = NKeyNU®) 

Since the measurements made from different local directions are independent, 
all measurements from 0's which fall into the bin can be easily combined to obtain 
compounded statistic of the measurement in the bin fi: 



m) _ Ee AUQ) - Ee KrM _ N on m - N b on (Q) QeQ 

(4.7) 

The critical value u c of the statistic U (Q) is set to five. If U(Q) is greater than five, 
the null hypothesis will be rejected and it will be said that the observed counts must 
have come from an astrophysical source. A measurement of the source strength can be 
performed. 

If the observed U (Q) is less than five, the null hypothesis will not be rejected and 
an upper limit corresponding to 2.3% error of the second kind will be made. A mea- 
surement of the source strength can be performed only if it is known (from other exper- 
iments) that the source exists. 

It should be remembered for probability interpretation according to the table 4.1 to 
be valid, the inequality (4.6) needs to be satisfied and for u c = 5 and typical value of 
a — 1/15 the number of events observed in the observation bin N on (il) should be about 
2- 10 6 . 
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4.6 Gamma-Ray flux measurement. 



Given the detector response to particles of different types and assumed source features, 
it is possible to predict the number of events N(Q) to be observed in the bin Q due 
to the source. Then, N(Q) can be compared with the actually observed number N(Q) 
and some statement regarding the assumed source features can be made. Indeed, let 
P(k, k, E, E, 0, 6, f, R) be the probability that the detector registers a particle of type 
k coming from local direction with energy E distance f from the apparatus and re- 
construction output information about the particle is k, 0, E, R. Then, the total number 
of events due to particles of type k to be observed from a region is: 

7V fc (fi)=V / P(k,k,E,E,e,e,f,R)F(k,E,e)T(Q)dEdEdededfdR 
r Jeen 



k 



where F(k, E, 0) is the number of particles of type k with energy E emitted by the 
source in local direction per unit area per unit time, T(0) is the time during which 
the source is located in local direction 0. The integration is performed over all possible 
values of energies E and E, all possible distances f and R and all directions in the field 
of view 0, but e Vt. 

The integration over core distances f, R, measured energy E and the summation 
over identified particle type k can be performed: 

P(k, E,Q,Q) = J2j P(k, k, E, E, 0, 0, f, R) dE dfdR 



N k (tt)= f F(k,E,Q)T(Q)\ f_ P(k,E,Q,Q)dQ 
J Jeen 



dEdQ 



loen 

If it is believed that the error in measuring event's direction does not depend on 
particle energy, P(k,E,Q,Q) can be factored as: 

P(k,E,Q,Q) = A k (E,Q)R k (Q\Q) 

The function A k (E, 0) is known as the "effective area" introduced in section 3.3.1; 
R k (Q\Q) is known as the angular resolution (or point spread) function. Then, the num- 
ber of events to be detected from the directions in the bin is: 



N k (n)= [ F(k,E,e)T(e)A k (E,e) \ f „R k (e\e)de 



dEdQ (4.8) 
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The integration should be performed over the exposure time to the whole source and 
given A k (E, Q)R k (Q\Q) can be done during data processing. Thus, by counting the 
number of events in an observation bin N k (Q) and comparing it with N k (Q), it is possi- 
ble to deduce some properties of the source function F(k, E, 6). If the null hypothesis 
is rejected, the difference N^Cl) = (N on (tl) — N% n (£l)) should be interpreted as 7-ray 
count. 

For instance, if a point source is considered with the source function F(7, E, ©) = 
F 5(E — E )5(Q — 6o(0) where 6 (i) is the source path in the local coordinates, the 
equation (4.8) gives: 



N s on (Cl) =Fj 5(0 - Q (t))T(Q)A,(E , 9) 17 . i? 7 (0|0) dQ 



dQ 



F ( A 7 (E , e (t)) \ I - R y (Q\Qo(t)) dQ 



dt 



F 



dt 



If the null hypothesis is not rejected, an upper limit corresponding to the error ( u can 
be set as (see equation (4.7)): 



N 



^OUtK^) 



leading to: 



F < 



«l(Cu)^ 


/Ee^(e) + Se^fy 


jA,(E ,e (t)) 







4.7 Optimal bin. 

If the Milagro detector had perfect angular resolution, then processing events from an 
infinitesimally small region of the sky around a point source would yield the maximum 
achievable sensitivity for source search as described in section 4.5. However, due to 
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detector's finite angular resolution, source events should be expected and accepted from 
some finite region around it. This, on the other hand, will increase the number of cosmic 
ray background events collected. Clearly, the optimal source acceptance region (called 
the "optimal bin") should have a configuration which provides the maximum power for 
the source search algorithm described in section 4.5. That is, for a given detector angular 
resolution function and given background event distribution on the sky, the optimal bin 
will provide the maximum value of the statistic U. 

In fact, the procedure for optimal bin construction, like the procedure for the con- 
struction of the best critical region, should be a part of the statistical test formulation 
and can not be modified based on observed data. Since, in the case considered here, the 
critical region on the values of the statistic U has been defined, the optimal bin construc- 
tion should be considered within the same framework. The optimal configuration Cl opt 
should maximize the value of statistic U(Q), thus the equation on the optimal search 
region Vt opt is: 

Sufi) _ Q 

Using the definition of U(Cl) from the equation (4.7): 



N b on (m + (E xe n %% M )/^n(^)] 



U{Q) = 

and neglecting ^-dependence of the term in the square brackets one arrives at: 

2 SN'Jfl) 1 5N b on (n) 



(4.9) 



The solution of this equation will be performed under the assumptions that the shape 
of the observation bin is circular with opening angle to « 1; that the bin is centered 
on a source occupying not more than a solid angle with opening o> << 1 and that the 
number of background events in the bin is proportional to its area (this is reasonable in 
the small angle limit). Then: 

NtM-^ and 1 2 



48 



true 
direction 



reconstructed 
direction 

/ 



source 




observation 
(p,()))\^ bin 



/ x 



Figure 4.4: Conceptual diagram of a small angular reconstruction error parameteriza- 
tion. 



Also, it will be assumed that the detector's angular resolution is described by a 2-D 
Gaussian with dispersion a 2 . This means that the error e on the reconstructed angle can 
be parametrized by two variables x and y which form a Cartesian coordinate system (see 
figure 4.4) and e = ^x 2 + y 2 . The probability of observing event (x, y) away from the 
true direction is: 



dR(x,y) = 



2na 2 



-(x 2 +y 2 )/2a i 



dxdy 



This representation is valid for small errors since the sphere may be treated as a plane. 

Under these assumptions, a point source with the source function 
F^(E, 6) = F(E)6(Q - 6 ) will produce the signal (see equation (4.8)): 



N s on = (t(0 o ) J E F(E)A 1 (E,Q )dE^ x £ 



-(x 2 +y 2 )/2<r< 



lx 2 + y 2 <u 2 2na 2 
= (T(Q ) j E F(E)A J (E,Q )dE^j x £± e ^ 2 ^ 2 de = 
= (V(9 ) J E F(E)A J (E : e )dE 



dxdy 



-0? /2a 2 



Then, the equation (4.9) becomes the equation on the optimal bin opening u>: 



{¥)' 



= 2 In 



1 + 



u opt ~ 1.585a" 



For a source whose source function F^(E, 6) is smooth within some opening angle 
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uj « 1 and zero outside and which is located in local direction O , the number of 
events detected from the source in the observation bin is (from equation (4.8)): 



A^(0 O ,^) = _ 

' x 2 +y 2 <u> 2 



x 2 +y 2 <u) 2 



T(0 O + e)F 7 (E, O + e)A(E, O + e) x 



x 



-({x-x) 2 +(y-y?)/2^ 



dxdy 



dxdydE 



where e = (x, y) — describes coordinate of a point inside the source bin relative to 
its center O . Introducing the notations uo a = uo/a and uj a = Co /a and substituting 
the coordinate system parameterization from Cartesian to polar as x = pa cos 0, y = 
pa sin the integration over can be performed and one arrives at: 



7V o s n (0 o ,u;W)=<7 2 



27r /-a) 



JO 



T(0 O + e)F 7 (£, O + e)A(E, O + e)dEx 



x 



pe 



-p 2 /2 



e-?/ 2 I (pp)pdp 



dp d(f) 



(4.10) 



where Iq(p) is the zeroth order modified Bessel function of the first kind and e = 

(pa cos 0, pa sin 0). 

In the expression (4.10) the integration over the angle can be performed by ex- 
panding the T(0 O + e)F y (E, O + e)A(E, O + e) in to Taylor series: 

T(0 O + e)F 7 (£, 0o + e)4(£, O + e) = T(0 o )F 7 (£, O )A(£, O ) + 

+e • V (T(0o)F 7 (E, O )A(0 O )) + 0(e 2 ) 

All even order corrections in e will give zero contribution to the integral (4.10) be- 
cause the integration is performed over 2n range. Keeping only the first order correc- 
tion in the source size uj a the number of events in the observation bin factors as: 



AP n (0 o ,cDW) = [a 1 J T(Q )F y (E,Q )A(E,e )dE) x 

x jT' pe^ 2 / 2 ( ^ e-?/ 2 I»{pp)pdp^ dp 
Hence, in the equation (4.9) the o -dependent factor will cancel and the optimal bin 
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Figure 4.5: Optimal bin size uo" t as a function of the size source co cr . 
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Table 4.2: Source size u a and corresponding optimal bin size u) a 



size will loose its dependence of the source location. The figure 4.5 shows the solution 
cu° pt of the optimal bin size problem (4.9) for the smooth source of size uj° ' . In the 
limiting case of the zero source size, the solution converges to the previously obtained 
Cj° pt = 1.585. It is also interesting to note that the fraction f° pt (Cu a ) of the signal events 
retained in the optimal bin is a very weak function of the source size: 

/°p'(0.0) = 0.715 /7*(1.05) = 0.721 /° p '(1.9) = 0.751 

The assumption that the number of background events in a bin is proportional to 
the bin's area is a good one for the Milagro data and the considered examples provide 
a good guideline for the choice of the observation bin size. Note that the constructed 
optimal bin is the optimal among the source centered circular bins and some other bin 
shape could be better. Nevertheless, the circular bin will be used in this analysis. 
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Chapter 5 



Photon flux at the Earth due to near 
solar neutralino annihilations 

'Yes,' I says to 'er, 'That's all very well,' I 
says. 'But if you'd been in my place you'd 
of done the same as what I done. It's easy to 
criticize,' I says, 'but you ain't got the same 
problems as what I got.' 
'Ah,' said the other, 'that's jest it. That's jest 
where it is.' 

George Orwell, "1984" 

While some introduction to the main goal of the current work has been done, the 
proper formulation of the problem is long overdue. 

There is overwhelming evidence that the Universe, and the galaxies in particular, 
are full of the "dark matter". There is no reason to assume that the Milky Way is any 
different. In this work, it is supposed that the weakly interacting particle (neutralino), 
predicted by super-symmetric theories, is the solution of the "dark matter" problem. If 
this is indeed the case, the neutralinos will form a halo around the Milky Way Galaxy and 
at the location of the Solar System the density of the halo neutralinos is often assumed 
tobep = 0.3 (GeV/cm 3 ). ' 

The neutralinos entering the Solar system may loose energy via elastic scattering 
with ordinary matter scatterers and become trapped in the Solar system. For simplicity, 
the Solar system is assumed to consist of the uniform density Sun only. This means that 
only the particles whose orbits cross the Sun can be captured on near-solar bound orbits 
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and their orbits will always cross the Sun. Due to the capture and repeated scatterings 
in the Sun, there will be a near-solar enhancement in the density of the neutralinos. 
This process is responsible for dark matter accretion in the Solar system. The dark 
matter diminution is due to neutralino-neutralino annihilations. The annihilation can 
not happen faster than accretion, otherwise, all dark matter would have annihilated by 
now. On the other hand, if accretion happens faster than annihilation, the Sun would 
constantly increase its mass. Thus, it is reasonable to assume that the Solar system 
has reached dynamic equilibrium and that the capture rate is equal to the annihilation 
one. Also, it is reasonable to assume that all possible elliptical orbits, not crossing 
the Sun have annihilated by now, as there is no mechanism to populate them by the 
incoming particles. Given that one of the annihilation channels is 7-ray production, 
one might expect an enhanced 7-ray signal from the neighborhood of the Sun due to 
neutralino annihilations. The purpose of this chapter is to estimate the 7-ray flux due to 
this process. 



5.1 General formulation of the problem. 

The problem of finding the density distribution of the particles in the Solar system can 
be addressed by kinetic theory. 1 If g(p, x) is the density of particles in phase space, then 
it should satisfy the Boltzmann equation 2 : 

>.) = ^* + ^4 = Cb fc.)] 

where C[g(p, x)] is the collision integral and the explicit dependence of g(p, x) on 
time has been dropped since in the considered model the process is assumed to be sta- 
tionary: dg(p, x)/dt = 0. The spatial density of particles n(x) is: 

n(x) = J g(p, x)dp 

The collision integral consists of two terms: one is due to neutralino annihilations 
C a [g(p, x)] and the other is due to scattering in the Sun C s [g(p, x)]. 

'The hydro-dynamic approach is not justifiable since the neutralinos do not interact with each other. 
2 Summation over repeated indices is assumed. 
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C S [#(P, a;)] = #0) / {^(P + 9, + q,x)- W(p, q)g(p, x)} dq 



1, \x\ < Rq 
0, \x\ > R Q 



e(x,R Q ) = 



where q is the particle momentum change in a collision, W(p, q) is the probability 
that a particle will change its momentum from p to (p — q) in a collision. 

A simplification can be made by noting that neutralino mass is much greater than that 
of any scatterer in the Sun and relative energy loss and momentum change of neutralino 
in a scattering are small. Thus, W{jp,q) is a quickly decreasing function of \q\ and 
diffusion approximation can be made: 



W(p + q)g{p + q,x) « 

d Id 2 
W(p,q)g(p,x) + qi—(W(p,q)g(p,x)^ + -qiqj j (W(p,q)g(p,x)) + 



±{A i (p)g(p : x) + -^ J 



d d 
C s [g(p,x)\ « — {Ai(p)g(p,x) + — (B ij {p)g(p,x))] 



where 



Mp) = J qiW{p,q)dq B^p) = ^J q i q j W(p,q)dq 



The function W(p, q) can be constructed by considering a structure-less elastic scatter- 
ing process where the angle of deflection of the incident particle in the center of mass 
reference frame is uniformly distributed between and it. 

The boundary condition for the problem can be formulated by assuming a Maxwellian 
distribution of velocities of galactic neutralinos. Then, in the Sun's reference frame the 
distribution will be shifted by the velocity of the Sun V in the Galactic disk: 

/ \ 3 / 2 (p-^xH)) 2 

lim g(p, x) = #00 (p) = Po - — 2~^~ e 2u o m * d 3 p 
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where m x is the neutralino mass. It will also be assumed that y2v$ = sJVq 
220 (km/s). 

The sought for annihilation rate density at a point x is simply: 



Needless to say, the task of solving the stated problem analytically or numerically 
is daunting even when W(p, q) has a simple structure. Direct computer simulations of 
the system will require enormous amounts of computer time. However, the distribution 
function g(p,x) is not the immediate goal of the project and only the distribution of 
the annihilation points is of interest. Therefore, it is proposed to perform computer 
simulations of annihilating particles only. This poses two problems: how to know that 
the particle will annihilate and what the boundary condition on the annihilated particles 
is. 



The stated problem is solved with "backward in time" simulation. The particles gen- 
erated at the annihilation points (so, it is known that the particles annihilated) are then 
propagated backward in time gaining energy in each scattering in the Sun until they 
exit the Solar system. If the distribution of the annihilation points is correct, the correct 
distribution of the annihilating particles at the boundary will be restored automatically. 
Thus, an algorithm should be devised to adjust an a priori distribution of the annihilation 
points in such a way as to reconstruct the correct distribution of annihilating particles 
at the Solar system boundary. This is possible if the particles carry information about 
their origin right until the annihilation point. The last assumption is reasonable since the 
neutralinos interact in the Sun only and their momentum does not change much in each 
scattering. 

Since two particles are required in the act of annihilation, the particle pairs will be 
considered. Let X describe the state of a pair of particles, then X will denote the 
pair state at the annihilation point and — the state at the boundary. The trajectory 
pairs can be divided into two classes: the first class containing the trajectories with both 
particles trapped in the Solar system via scattering in the Sun and the second — the 
trajectories when at least one of the trajectories was not trapped and annihilated in its 
first flight through the Solar system. Supposing that the time of flight through the Solar 





5.2 General idea of the solution. 
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system is much smaller than the mean life time of the particle, the contribution from the 
second class can be neglected and only captured trajectory pairs can be considered. Let 
the joint probability that a pair of captured particles annihilated at X and entered the 
Solar system at be P(X OQ , cap, X ). The task of the project is to find the distribution 
of the annihilation points or P(X , cap) and according to Bayes' theorem, the joint 
probability P(X 00 , cap, X ) is: 

P(Xoo, cap, X ) = P(X \X oo , cap)P(X OD , cap) = P{X Xl \cap, X )P(X , cap) (5.1) 

J P{X O0 , cap, X ) dX = P(X oc , cap) = J P^X^cap, X )P(X , cap) dX (5.2) 
and 

P(X 00 ,cap) = P(cap|X 00 )P(X 00 ) 

where P(X 00 ) is the known distribution of all trajectories entering the Solar system 
at infinity and P(cap|Xoo) is the conditional probability that the trajectories X^ will be 
captured. 

The probabilities P(X 00 ) and P(cap|Xoo) can be found analytically (see 
sections 5.2.3 and 5.2.4 correspondingly) and P(X 00 \cap, X ) can be constructed using 
a backward in time computer simulation (see section 5.2.2). Also a method of solving 
the equation (5.2) for P(X , cap) is described in section 5.2.1. 



5.2.1 Solution of the Fredholm equation. 

The equation (5.2) is the Fredholm equation of the first kind with respect to the sought 
for function P(X , cap): 

P(X 00 ,cap) = J P(X oo \cap,X )P(X ,cap)d(X ,cap) (5.3) 

The kernel P(X OQ \cap, X ) of this equation is not known analytically and is con- 
structed in simulations (see section 5.2.2). It should be noted that due to the nature 
of the kernel construction method, the kernel P(X OQ \cap, X ) represents a "list" of 
"(X , cap)"s and "(X^, cap)"s which are "connected" by the simulation process. Be- 
cause the initial state (X , cap) is random and the propagation process is stochastic, it is 
statistically improbable to have repeated pairs in the list. This implies that all knowledge 
regarding the kernel can be expressed as: 
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P(X oo \cap,X ) 



1 , { (Xoo , cap) , (X , cap) } e list 
0, {(Xoo, cap), (X , cap)} £ list 



Thus, the equation (5.3) becomes: 

P(X Q , cap) = P(Xoo, cap) 

In other words, when annihilation state (X , cap) of a captured particle is associated 
with a state at the boundary (X^, cap) by the propagation process from section 5.2.2 
the relative contribution P(X , cap) from the state (X , cap) is P(X oc , cap). 

Thus, generating the initial states (X , cap) such that all final states (X^, cap) are 
sampled will lead to the solution of the equation (5.3). 

To construct the histogram of the radial distribution of the annihilation points it 
should be noted that the state X = (r , v 1: v 2 ) is the position and the velocities of the 
two particles at the annihilation point, thus: 

P(f ,cap)— / P(X ,cap)dv 1 dv 2 = / P(f , v 1 , v 2 ) dv ± dv 2 

Joj(fo) Jw{?o) 

where Lu(r ) is the velocity volume over which the integration is being performed. This 
volume should include all particles which are captured and is finite. If the above inte- 
gration is performed by the means of the Monte Carlo method with uniform sampling 
in the velocity volume P(r , cap) becomes: 

P(f ,cap) = U ^\ J2 P(r ,vi,v 2 ) 

where N v (f ) is the number of sampled points. 

Because the histogram of the annihilation points H(f ) is defined is the average of 
P(f , cap) in the f bin, H(f ) is: 

* (ft) - £ E nf a ,ca P) = -L ? g| e P (r -„, a) 

where N ro is the number of entries in the r bin. 

The obtained expression for the histogram H(f ) of the annihilation points may be 
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simplified by choosing a fixed large volume of the velocity space co(f ) = Q and by 
noting that in any finite random sample the probability to observe two different pairs of 
orbits passing through the same point f is zero leading to N v (f ) = 1. 

5.2.2 Transition tables { (X^, cap) , (X , cap) }. 

The transitions from to X for captured particles are found with backward in time 
simulations and was mentioned before. Because the particles evolve independently, each 
particle is propagated from its initial state (annihilation point) backward in time until it 
encounters an interaction in the Sun. At this point, momentum is changed according 
the rules of elastic scattering, the energy is gained, and the new angular momentum is 
computed. Afterwards, the particle is propagated until it encounters the next scattering. 
This process repeats until the accumulated energy is greater than zero, which means that 
the particle is no longer bound to the Solar system and its state at infinity is found and 
recorded. 

Because the particle may spend long periods of time between interactions in the Sun, 
one must solve the equations of particle motion analytically and use the results. The 
motion in central potential fields is integrable and each trajectory is defined by integrals 
of motion: the total energy of neutralino E and its angular momentum J. However, the 
kinematics depends only on mass densities of these quantities: S = E/mand J = J /m 
and mostly those will be used in the calculations. 

The major simplification in the trajectory calculation comes from the fact that the 
energy and angular momentum are conserved between the scatterings. This means that 
between the scatterings the motion is executed in one plane and only rotation of the 
whole orbit is possible. The Runge-Lenz vector JC, which fixes the orientation of the 
particle orbit and is conserved outside the Sun may change its direction only when par- 
ticle passes through the Sun. Since between the scatterings the trajectory inside the Sun 
does not change its properties, the particle passage through the Sun can be tracked by 
rotation of the Runge-Lenz vector. Any point on the particle trajectory can be specified 
by its angle relative to the current direction of the Runge-Lenz vector. 

The trajectory length inside the Sun plays the main role in the propagation process 
because it is the quantity which defines when the next scattering should occur. Given 
an initial point and the path-length inside the Sun until the next scattering, the point of 
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next scattering can be found by rotating the Runge-Lenz vector on the appropriate angle. 
Thus, the task is to convert the trajectory length inside the Sun into the angle of rotation. 
This problem can be solved analytically for the specified solar model. 

The act of scattering can also be described as rescaling and the rotation of the ve- 
locity. Thus the whole process of particle propagation from its annihilation point back 
to infinity can be described as a sequence of rotations applied to the Runge-Lenz vector 
and the velocity rotation with rescaling. 

The details of the simulation process are described in the appendix D. 



5.2.3 Distribution of neutralinos at infinity P(X OQ ). 

Because the state describes the state of two identical particles at the boundary, the 
P(X OQ ) function is constructed as the product of two identical distributions describing 
a single particle: 

P(Xoo) = 5-00(^1)5-00(^2) 

The expression for g^x) is a simple generalization of the results obtained earlier 
in [24] and [43]. However, a self-consistent derivation is provided here for complete- 
ness. 

Let us choose a sphere of a large radius R around the Sun so that the effects of Sun's 
gravity are negligible and the velocity distribution of the particles is a known function 
f(v)d 3 v and does not depend on the point on the sphere. Let n x be the concentration of 
particles at the sphere. The number of particles entering the Solar system per unit time 
with velocity v from the surface element dA is then: 

dN = n x f{v){v-dA)d 3 vdt 

where we are interested in the particles with (v ■ dA) < since the particles should enter 
the sphere. 

_ — * — * 

Since we are considering the sphere dA = R 2 sin 9 d6 d(j)j and dA ] \ R we can 
choose the coordinate system (9, <fij) so that 9 is counted from the direction of the 
velocity v, then: 



[v ■ dA) = (v ■ R)R sin 9d9d(j)j = ^-d {v 2 R 2 sin 2 6) d(f) 



dj 2 d(pj 
2v ~ \ 2v 
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Then, the number of particles entering the Solar system is: 

dN = n x f(v) — {v-R)<0 

If we are interested in the total distribution, we must note that since the velocity 

— * 

distribution does not depend on the spatial point, for every particle with (v ■ R) < 

— * 

there will be exactly one particle with (v ■ R) > 0. Hence, the number of particles 
entering the Solar system with the velocity v and magnitude of angular momentum J 
per unit time is: 



xJK ' 4v ' v ; dtd 3 vdj 2 d<pj 4v v 

Consider the case when the velocity distribution at infinity is spherically symmetric 
as in [43], 

f(v)v 2 dv = An(2na 2 )- 3 / 2 e~ v2 ^ 2 v 2 dv 

then after integration over the spherical coordinates of the velocity and dtfij one arrives 
at: 

dN 
~dt 

This is the formula (2.7) from [43]. 



2n 2 n x f(v)v dv dj 2 = 2n 2 n x f(v) dS dj 2 



If the motion of the Sun with speed V in the locally isotropic Galactic halo is taken 
into account, then, as in [24]: 

f( v y d v = f f(v)— = s J^^ e -v V^ Un(2na 2 )- 3 / 2 e- v2 ^ 2 v 2 dv 

and the rate at which the particles enter the Solar system per angular momentum per 
speed is: 

r d 3 V 7T ~ 

dN = 7in x dtdJ 2 / f(v) = -n x dtdj 2 f(v)v 2 dv 

This is the expression which will be used for calculation of the distribution of the parti- 
cles at infinity. 
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5.2.4 Capture probability P(caj9|X QO ). 

Because the state describes the state of two identical particles at the boundary and 
because each particle is captured independently, the P^aplX^) function is constructed 
as the product of two identical capture functions describing a single particle: 

P(cap|Xoo) = g ca p(xi)g cap (x 2 ) 

Also it will be assumed that the capture happens in one collision. In other words, 
after the first collision (forward time) the particle will have negative energy 3 . This will 
greatly simplify calculations, while higher order corrections can be considered. The 
argument for this is that the mean-free-path of neutralinos in the Sun A = l/n p a PX is of 
the order of 10 4 — 10 9 (Rq) for expected values of a px . The energy loss in a collision is 
(see section D.9): 



AS = Wy™ <>) ± = 2*7(1- cos _ x = , _ 



(r] + l) 2 2 (ri + 1) 2 V V K '> ' m v 

where £ is the energy before the collision and U (r) is the potential energy at the colli- 
sion. Since 9 is the scattering angle in the center-of-mass system and no cross-section 
structure is assumed, cos 9 is distributed uniformly between —1 and 1, this leads to the 
fact that v is uniformly distributed between zero and j^p^p. ■ 

( v > c , need for capture 

Safter = £-A£<0 £-v (£-U (r ))< { . £ - u{r \ v „ , 
a}ter v v 11 \ < v < 7-^2 allowed range 

Thus, the probability that a particle will be captured in one collision at the distance r 
from the center of the Sun is: 



4r? 



4r? f 



(r? + l) 2 £-U{r) 



e 



4?7 5 



(?7 + l) 2 S-U(r) 




The probability that a particle will travel distance y without scattering and scatter in 



3 The particle may scatter twice or more in the first pass through the Sun, but it is assumed to become 
trapped after the first collision. 



61 



y,y + dy is: 

dp = \e~y/ x dy » % 
A A 

Again, this approximation is valid since the mean free path of neutralinos inside the Sun 
A is much greater than the particle trajectory length inside the Sun. 

The probability that the particle with energy £ and angular momentum J will loose 
energy in one collision to become captured on the bound orbit is: 

fL(£,J) 1 

g C a P {x) = jPca P (£, r{y))dy 

where L(£, J) is the path-length inside the Sun. 
From the energy conservation law: 

v 2 r 2 J 2 

§ = V '2[£-C/(r)], £ = ^2\£ - U(r)] - J*/r* 



Qcap{%) 



2 f^max 



2[£ - U{r)\ 



\Jr min \2[£-U{r)]-J 2 /i 



Arj 



Arj 



£ 



(77 + 1) 2 £-U(r) 



dr 



The r min is the minimal distance from the Sun's center to the orbit and can be found 
from the equation of motion: 

£ = J 2 /Zr 2 mm + U(r mm 



72 /2r 2 

' ""' """ 2R . 



© 



(3 - r 2 mm /Rl 



r ■ = R z 

' mm * l O 



(3aR Q + 2£R 2 Q ) - J(3aR Q + 2£R%) 2 - AaR Q J 



2aR, 







The r max is computed from the restriction on maximum capturable energy and 
should not be greater than R Q or less than r min . 



A-q 



(77 + 1) 2 £-U{r) 



> 



£ < 



A-q 



(77-1) 



Mr) 
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r 2 -P*U ^- l ) 2 c\ r . < r <R 5 c 4 ^ 3a 

W - ^0 ^ a 2r? ^ , rw < w < it; , fc < ^ _ 1)2 2i?£ 







The integrals which need to be executed to find g cap (x) are elliptical integrals. g cap (x) 
can be written in the form: 



1 / a-6y 

W^) = t / 9 \ — — j— 2 — ~dy - 



A Jr 2 . V ay — by 2 — c 
(77 + l) 2 8 r 2 ™™ dy 



t r 

A Jr 2 



2 V ^Jr 2 mm ^(ay-by 2 -c)(a-by) 

where 



a = 28 + —, b= w , c = J 2 



From Gradshtein and Rizhik(Russian 3.141-2 page 245): 



From Gradshtein and Rizhik(Russian 3.131-3 page 233): 



f^max j 1 1 f r max I j_ 

Xi.„ V (ay-by 2 -c)(a-by) dy= bJc ]j (y - A)(y - B)(y - C) dy 

EF(>y,q) 



lrL n Uay-by 2 -c){a-by) y b Jc ]/ (y - A)(y - B)(y - C) 

2 



bJ(A-C) 



• r 2 max -C B-C 

7 = arcsm V^^' g = V^ 

Further simplification comes from the fact that A — C = B 
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m x , (TeV) 


fout 


I(m x ) X 10 -43 cm 2 o .3GeV/cm 3 ' (" ) 


0.1 


0.195 


1.65 10 18 


0.2 


0.195 


4.17 10 17 


0.5 


0.196 


6.72 10 16 


1.0 


0.199 


1.68 10 16 


2.0 


0.201 


4.22 10 15 


5.0 


0.2 


6.72 10 14 


10.0 


0.2 


1.69 10 14 



Table 5.1: Summary of the simulation/computation results. Capture integral I and the 
fraction f out of annihilations between 1 and 2 solar radii as a function of neutralino mass 

m x . 

Where the elliptical integrals are: 

EFU,k)= 1/Jl -k 2 sin 2 (t)dt = / , |Jfe|<l 

J ° V Jo ^/(l-t 2 )(l-k 2 t 2 ) 



EE((j),k) = J^l-k 2 sm 2 (t)dt = J i _ dt, \k\<l 

5.3 Predicted photon flux. 

The results of the computer calculation are summarized in the table 5.1 for several se- 
lected neutralino masses. About 40 — 50% of particles annihilate outside the Sun, but 
their distribution is a sharply falling function of distance from the Sun (see figures 5.1 
and 5.2), so only the annihilations happening between one and two solar radii will be 
considered to produce detectable signal. The fraction of neutralinos annihilating be- 
tween one and two radii of the Sun is denoted as f out in the table. 

Only a small fraction of the annihilated particles will produce photon signal. If pho- 
ton yield for producing a photon with energy E 1 per neutralino in neutralino-neutralino 
annihilation is b 1 (E 1) m x ), the total number of photons produced per second is: 

rf$o = H m x) ■ fout ■ & 7 (^ 7 , m x ) dE 1 
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Some of the produced photons will be absorbed by the Sun. The fraction of the 
photons escaping the Sun is f eS ca P e and is of the order of 1/2. The distance between the 
Earth and the Sun is L e which leads to the flux of number of photons per unit area per 
time at Earth from neutralino annihilations as: 

dF x( E l) = H m x) ' font ■ by( E li m x) ■ f escape/ dE 1 



dF x (E 7 ) = 

PO &px 1/771 \ f° ut ' f escape I\ m x) j 771 / -2 —1\ 



0.3 {GeV/cmt) 10~ 43 cm 2 7V n x > 0.2-1/2 2.8 • 10 28 

(5.5) 

The photon yield may have the following structure: 

b,(E„m x ) = ^K)5(S 7 - m x ) + b;(m x )P(^-) 

where P{J^-} is the probability to produce photon with energy E 1 , E y < m x in the 
continuous spectrum neutralino to photon annihilation process. 

There are indications (see [9]) that the continuum spectrum probability has the form: 



^m x ' m x ^m x ' 
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Figure 5.1: Radial distribution of the annihilation points for m x = 200 (GeV) and 
a px = 10 -43 (cm 2 ). Vertical scale is in arbitrary units, horizontal scale is in R Q . above 
25 • 10 6 neutralino annihilations 
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Figure 5.2: Radial distribution of the annihilation points for m x = 1000 (GeV) and 
a px = 10 -43 (cm 2 ). Vertical scale is in arbitrary units, horizontal scale is in R & . above 
22 ■ 10 6 neutralino annihilations 
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Chapter 6 



Outcome of the test for presence of the 
photon flux from the Sun and its 
implications 

His courage seemed suddenly to stiffen of its 
own accord. 

George Orwell, "1984" 

The gamma ray signal from neutralino annihilations near the Sun should appear as 
an excess number of events from the direction of geometrical center of the Sun over the 
expected background. Observation of the Solar region can be performed by tracking the 
Sun on the Celestial sphere using the one-arc-minute precision formulae for the Sun's 
Celestial coordinates from [12]. The interpretation of the observed signal, however, is 
not an easy problem. Largely, this is due to the fact that the cosmic ray background is 
not expected to be uniform; the Sun absorbs the cosmic rays impinging on it and forms a 
cosmic-ray shadow. The situation is complicated by the magnetic fields of the Earth and 
the Sun. Due to bending of charged particles trajectories in magnetic fields, the Sun's 
shadow will be smeared and shifted from the geometrical position of the Sun in the TeV 
range of particle energies. On the other hand, in the presence of strong Solar magnetic 
fields, lower energy particles can not reach the surface of the Sun and are reflected from 
it. Such particles are not being removed from the interplanetary medium and may not 
even form a cosmic-ray shadow of the Sun. Therefore, it is difficult to ascertain the 
exact shape of the cosmic-ray shadow at the Sun's position and deduce excess above it. 
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N 

ly Oil 


N b 

on 


U 


Sun 


137211 


137728 


-1.35 


Moon 


49762 


50064 


-1.31 



Table 6.1: Number of events in the optimal bin centered on the Sun and the Moon (see 
section 4.5). 

The effect of the Earth's magnetic field and the Solar wind can be studied by ob- 
serving the shadow of the Moon during solar day. If the solar magnetic field is weak, 
the shadows of the Sun and the Moon should be very similar because of the geometry 
of the problem. The Sun and the Moon cover similar size regions on the celestial sphere 
and traverse similar paths on the local sky in one year of observation. In addition, the 
Moon is far enough from the Earth to be considered outside of the effect of the Earth's 
magnetic field, so is the Sun. 

6.1 The data set. 

The data to be used in this work was chosen to satisfy the following conditions: online 
reconstruction, the number of photo-tubes required to trigger the detector greater than 
60, the number of photo-tubes used in the angular reconstruction (Nfit) greater than 
20, zenith angles smaller than 45 degrees and all events should pass the gamma/hadron 
separation cut. The data used were collected between July 19 2000 and September 10 
2001. The dates are motivated by introduction of the hadron separation parameter into 
the online reconstruction code on July 19, 2000 and detector turn-off for major repairs 
on the 11th of September 2001. Several data runs were disregarded from the dataset 
which included calibration runs and the data when the online DAQ was in an unstable 
regime. 

For the Sun analysis a ±5° regions around the Moon and the Crab nebula were 
vetoed from the data set. For the Moon analysis, same size regions around the Sun and 
the Crab were vetoed. Overall, 1 164.7 hours of exposure on the Sun and 423.5 hours of 
exposure on the Moon during the day time is obtained in this data set. The number of 
events in the optimal circular bin of 1.26° radius centered on the Moon and the Sun is 
given in table 6. 1 and the sky maps with corresponding exposure graphs are presented in 
figure 6.1. The sky maps are generated according to the equation (4.1) with the vertical 
axis pointing to the Geocentric Geomagnetic North dipole pole. 
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Figure 6.1: Significance maps of the regions of the sky around the daytime Moon(left) 
and the Sun(right) and the corresponding source exposure as function of zenith angle in 
hours per degree. The color code is the value of U (see equation (4.7)). 
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m x (TeV) 


A (cm 2 s) 


S (cm 2 s) 


0.1 


1.055- 10 11 


0.000 


0.2 


8.772 • 10 11 


4.969 • 10 7 


0.5 


6.070 • 10 12 


2.634 • 10 y 


1.0 


3.389 • 10 13 


2.127- 10 10 


2.0 


1.600 • 10 14 


1.280 • 10 11 


5.0 


5.942 • 10 14 


1.208 • 10 12 


10.0 


1.608 • 10 15 


5.575 • 10 12 


20.0 


3.684 • 10 15 


2.136- 10 13 


50.0 


8.030 • 10 15 


1.035 • 10 14 



Table 6.2: Coefficients of the flux limit calculation (see equation (6.1)). 



m x {TeV) 


F s < (cm-is- 1 ) 


a PX PO Ud ^ 
10~ 43 cm 2 0.3GeV/cm 3 7 


0.1 


4.54- 10~ 8 


770 


0.2 


5.46 • 10~ 9 


351 


0.5 


7.89 • IO" 10 


328 


1.0 


1.41 • IO" 10 


234 


2.0 


2.99- IO" 11 


204 


5.0 


8.06- 10~ i2 


339 


10.0 


2.98- 10~ 12 


512 



Table 6.3: The upper limit on the monochromatic photon flux due to near- solar neu- 
tralino annihilations and corresponding upper limit on the o- px p 6^. 

6.2 A limit on possible gamma-ray flux due to near-Solar 
neutralino annihilations. 

Because the shape of the solar shadow is not known, not to claim a false signal the null 
hypothesis is formulated as cosmic-ray background is uniform and there is no 7-ray 
emission from the solar region. The mean value of the statistic U (see equation (4.5)) 
is equal to zero under this hypothesis. Based on the results of the measurement (see 
table 6.1) the formulated null hypothesis can not be rejected with significance 2.867 • 
10~ 7 (see table 4.1) and an upper limit on the possible 7-ray flux from the solar region 
should be obtained. 
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Figure 6.2: The values of (F$,F C ) below the lines are allowed based the constructed 
upper limit for corresponding neutralino masses. 



m x (TeV) 


F c < (cm~ 2 s x ) 


l -43 cm 2 .3GeV/cm :i u l 


0.1 






0.2 


9.64- 10~ 5 


6.47- 10 6 


0.5 


1.82 • 10~ 6 


7.58 • 10 5 


1.0 


2.25 • 10~ 7 


3.75 • 10 5 


2.0 


3.74- 10- 8 


2.48 • 10 5 


5.0 


3.97- 10~ 9 


1.65 • 10 5 


10.0 


8.59- 10- 10 


1.42 • 10 5 



Table 6.4: The upper limit on the continuum photon flux due to near-solar neutralino 
annihilations and corresponding upper limit on the a PX p b^. 
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The deficit of events from the direction of the Sun can not be greater than that pro- 
duced by the Moon because of Sun/Moon similarities. To be conservative in setting 
the upper limit, the strongest event deficit produced by the Moon in 5° radius from its 
position should be used as a correction for possible presence of the shadow of the Sun. 
The smallest value of the statistic U observed in the sky map centered on the geomet- 
rical position of the Moon is —3.3 (see figure 6.1). The exposure on the Sun is about 
2.75 greater than that on the Moon, leading to the estimated maximal deficit in the Sun's 
direction computed in terms of U as: U™™ = —3. 3\/2. 75 = —5.5. 

Thus, the upper limit on the photon flux from the region of the Sun corresponding to 
the significance 2.867 ■ 1CT 7 with error of the second kind 2.275 • 1CT 2 (see table 4.1) is 
computed based on the value of the statistic U : 

Ul = 5.0 + 5.5 + 2.0 = 12.5 

leading to the upper limit on the mean number of the gamma counts: 

N <N U = Ul ^N b + aN s = 4791 

The differential photon flux due to neutralino annihilations taken from [9] has the 
form of: 

dF(E) pA(p r.(£ > 0.01) (JLy 3 ' 2 e-^ 

where Fs is the integral flux due to a 5-function-like photon annihilation channel and 
Fc{-^~ > 0.01) is the integral flux for > 0.01 due to continuum photon spectrum 
annihilation channel of neutralino s with mass m x . 1 

Computing the number of events to be observed by the detector using the for- 
mula (4.8) for the given spectrum and following the procedure for setting an upper limit 
(from section 4.6), it is possible to set a constraint on the integral fluxes in the form: 

F s - A + F c -S < 4791 (6.1) 

Fs and F c are in the units of cm~ 2 s~ l and the coefficients A and E are given in the 
table 6.2 for different neutralino masses. The figure 6.2 shows the region of parameter 



1 J x -3/2 e -a Xdx = _2e_^_ _ 2 ^Erf(y^) 
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Figure 6.3: The values of yp a px b^, p crp X b^J below the lines are allowed based the 
constructed upper limit for corresponding neutralino masses. 

space (F c , F$) restricted by the upper limit. 

The upper limit on the monochromatic photon flux due to neutralino annihilations 
and corresponding limit on <7 px p &^ (see equation (5.5)) are presented in table 6.3. The 
upper limit on the continuous photon flux with energies above 0.01m x due to neutralino 
annihilations and corresponding limit on a PX p b^ (see equation (5.5)) are presented in 
table 6.4. The figure 6.3 depicts the combined limit on a PX p b^ and a PX p b^. 
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Chapter 7 
Conclusion 



The landscape that he was looking at re- 
curred so often in his dreams that he was 
never fully certain whether or not he had seen 
it in the real world. 

George Orwell "1984" 

The Milagro data set collected during 2000-2001 has been analyzed and searched for 
the evidence of a steady near-lTeV 7-ray flux from near-solar neutralino annihilations. 
As a result of the analysis, it was argued that no evidence for the gamma-ray signal due 
to such a process has been found. The upper limit on the gamma-ray flux due to such a 
process with significance 2.867- 1CT 7 and the power (1 —2.275- 1CT 3 ) has been set. Even 
in the absence of a clear signal the constructed upper limit may constrain the values of 
free parameters of supersymmetric models. 

The interpretation of the constructed limit on the gamma-ray flux is highly model 
dependent. It is based, for instance, on assumptions regarding the shape of the velocity 
distribution of the dark matter in the halo and the assumed structure of the Solar System. 
The current work presents a calculation of the neutralino annihilation rate density as a 
function of distance from the Sun and the neutralino capture rate onto near-solar bound 
orbits. It is shown that in the considered model only about 50% to 60% of annihila- 
tions happen inside the Sun. The calculation allowed translating the established limit 
on the gamma-ray flux from neutralino annihilations to the limit on the product of the 
neutralino-proton scattering crossection a px , the integrated photon yield per neutralino 
in neutralino-neutralino annihilation 6 7 and the local galactic halo dark matter density 
Po- 
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To the knowledge of the author the current work presents a first attempt to set a con- 
straint on the parameters of supersymmetric models by observing high energy gamma 
rays from the region of the Sun. Continuous improvements in reconstruction algorithms, 
detector modifications and longer observation times will led to a better upper limit. One 
of the factors which lead to a deterioration of the constructed upper limit is the inability 
to compensate for presence of the Solar cosmic-ray shadow due to the intricate struc- 
ture of the Solar magnetic fields. Once the cosmic-ray shadow of the Sun is understood 
quantitatively, it may be possible to improve upon the limit. 
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Appendix A 
Poisson distribution 



A.l Definition 



The Poisson distribution arises in many problems as the distribution of the number of 
occurrences of some event over an interval of time or region of space. The distribu- 
tion assumes that an event can occur at random at any time or point in space and that 
the probability of event occurrence does not depend on any other event. The Poisson 
distribution is defined as: 



P (*;A) = ^e-* 



(A.l) 



oo oo \ k 

£p(M) = e-*£-=e , 

k=0 k=0 K - 



<k>=Y: kp(k; A) = e- x ]T 



fc=0 



A;=0 



-A 



oo 



D(k) =< k 2 >-<k> 2 = e- A £ 



kJAk-l)\ 

oo k 2 X k 



A 



k=0 



k=i 



A 



(k- l)\ k 

(k - 1)! + (k - 1)! 



A 2 = e" 



k\ 



A 2 e A + \e ; 



X 2 = 



A 2 = A 
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A.2 Gaussian Limit of Poisson Distribution 



Substituting the k\ in the Poisson distribution (A.l) by the approximate expression using 
the Stirling formula: 



the Poisson distribution becomes: 

p(M) 



eA 



or 



In 



k In 



2ti k \ k 

eX 



k 



p(k;\)V2irk 
Let k = A + 5 than 

ln[p{k;X)V2^k] « (A + 5) 



X = k 



1-ln^ 
A 



A 



1 - In 



A_+5 
A 



-A 



(A + 5) 



Thus: 



p(fc;A) 



1 



2tt (A + 5) 



e 2A . eeiA 2 



For sufficiently small 5 such that | << 1 and 



<< 1 the Poisson distribu- 



tion approaches the Gaussian distribution with mean and dispersion equal to A: 



e^ 1 — 1 



« 1 =>• 



<5 3 



6A 2 



«1 -«\ -<l, VA>6 

A V A 



1 (fc-A) 2 

e 2 * 



2ttA 



(A;-A) J 



6A 2 



« 1, VA > 6 



(A.2) 
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Appendix B 
Calibration 



For, after all, how do we know that two and 
two make four? 

George Orwell, "1984" 

The desire to reconstruct the position of events on the Celestial Sphere with system- 
atic errors much less than 1° dictates that the times registered by PMTs have to have 
resolution about 1 (ns) and the locations of the photo-tubes be known to about 10 (cm) 
accuracy in horizontal and about 3 (cm) in vertical directions. To meet the latter require- 
ment photographic and theodolite surveys of the pond were performed. At the end of 
the construction period, when the pond was filled with water, an "as-built" measurement 
of the PMT elevation was done. 

Even though great care has taken to construct all PMT channels as uniformly as 
possible, remaining systematic differences in PMT channels should be studied and re- 
moved. This includes synchronization of all TDC modules (find TDC time offsets and 
conversion factors) and compensation for the PMT-pulse amplitude dependence of TDC 
measurements (known as the slewing correction). 

To correctly reconstruct the shower front, shower size and, ultimately, to estimate the 
energy of the primary particle, the relative "pulse-height" to photo-electron (PE) con- 
version must be determined to interpret all PMT amplitude measurements in a common 
unit for each event. This is then translated to an absolute scale of the energy deposited 
in the water. 

Full description of the calibration system with its components, operation, data anal- 
ysis and other comments is available in [13] and references cited therein. 
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Figure B.l: Calibration system setup 



B.l Calibration system setup. 

The Milagro calibration system has been designed to reflect all the above goals and is 
based on the laser - fiber-optic - diffusing ball concept used in other water-Cherenkov 
detectors (See, for instance, [6]). A computer operated motion controller (Newport 
MM3000) drives a neutral density filter wheel to attenuate a pulsed nitrogen dye laser 
(Laser Photonics LN120C) beam. The beam is directed to one of the thirty diffusing 
laser balls through the fiber-optic switch (DiCon MC606) as shown on Fig B.l. Part of 
the laser beam is sent to a photo-diode. When triggered by the photo-diode, the pulse- 
delay generator (Stanford Research DG535) sends a trigger pulse to the data acquisition 
system. A laser fire command is issued by the motion controller, providing full automa- 
tion of the calibration process. The balls are floating in the pond so that each PMT 
can register signals from more than one light source. Such a redundant setup allows 
calibration of the PMTs and the electronics. 

Calibration data was collected by stepping the filter wheel through a full circle in 
10 degrees increments for each laser ball. On each laser ball - filter wheel setting about 
2000 laser triggers at 20 (Hz) and about 1600 "random" triggers (with no light input) 
at the rate of 400 (Hz) were sent to the data acquisition system. Raw data from all 
PMT channels was recored and analyzed. Only 2- and 4-edge events with correct po- 
larity were selected to ensure proper ToT determination. For purposes of occupancy 
measurements all data was used without any edge selection. 
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B.2 Timing calibration. 



The importance of accurate time readings from the PMT channels can not be over 
stressed as the quality of event reconstruction depends on it. The issues which need 
to be addressed are described in this section. 

B.2.1 TDC Conversion Factor. 

The time of PMT pulse threshold crossings is read by LeCroy 1887 FASTBUS TDC 
modules. These digital devices measure time in the units of "counts" and, according 
to specifications, each count corresponds to 0.5 nanoseconds. Introduction of known 
variable delays in the calibration-DAQ trigger logic 1 allowed observation of common 
time shifts in all PMT channels to verify the TDC conversion factors at 2.0000 ± 0.0003 
counts per nanosecond. This assured that all TDC modules operate on a common scale 
and allowed for a simple interpretation of TDC measurements. 

B.2.2 Electronic slewing correction. 

Time response of a PMT channel depends on the input light intensity and is called 
electronic slewing. Indeed, a weaker pulse will cross the discriminator threshold later 
than a strong one arriving at the same time (see figure B.2). Based on the ToT PMT 
pulse model described in 3.1.2, when the size of the pulse is described by the time 
over threshold, a slewing correction can be devised by studying the PMT pulse arrival 
time (T start ) as a function of ToT for different filter wheel transparency settings. The 
time of light pulse emission is supplied by the photo-diode and is believed to be free 
of slewing effects since the light level incident on the diode is constant. The slewing 
correcting curve is found by fitting the obtained calibration data with a polynomial (see 
figure B.3). 

Note, however, that obtained T start includes time propagation of the laser pulse in 
the detector medium and optical fibers. These may vary from PMT to PMT and should 
be taken into account to produce true T start versus ToT dependence. 

The procedure described above should be performed independently for both thresh- 
old levels HiToT and LoToT . 



Special TDC calibration data should be taken for this procedure. 
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Tstart 



[Tstart 




Figure B .2: Illustration of the electronic slewing. Stronger pulses cross the discriminator 
threshold earlier than the weaker ones. 




Figure B.3: Plots show T start vs ToT data obtained for calibration (left) and the poly- 
nomial fit to the data (right). The units of both axes are TDC counts. 
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B.2.3 Speed of light in water, fiber delay. 



In order to correct for propagation time of light in the detector, coordinates of PMTs 
and laser balls as well as the speed of light in water must be known. PMT and laser ball 
coordinates are known from the survey. Only a typical index of refraction of water is 
found in reference tables and fiber optic delays may vary from laser ball to laser ball, 
thus, all these parameters need to be measured with the calibration system. 

Interestingly enough, this problem can be easily solved [17] if several PMTs can 
register light from two laser balls (cross-calibration). The times measured from two 
different laser balls T} tart and T^ tart on the same PMT after slewing correction should be 
identical. The non-zero difference between the times T] tart and T^ tart can be attributed 
to an error in the water propagation time ^.propagation and/or difference in the laser balls' 
fiber optic delay A/j ber . If t is defined as: 



' = 1 start ~~ 1 start ~ ^ fiber ~ ^propagation 

it will be zero in absence of errors in water propagation time and fiber delays. 

The distribution over all PMTS of observed r from a given laser ball pair can be 
constructed and studied. Since the A f iber is constant for the given laser ball pair and 
^■propagation depends on the relative PMT ball positions, the use of correct speed of 
light will yield the minimal width of the r distribution. After that, the mean of the 
distribution can be interpreted as the fiber optic difference Af iber . Note that a PMT in 
close proximity to any one of the laser balls will have enhanced sensitivity to speed of 
light variation, while PMTs located half way between the laser balls will have enhanced 
sensitivity to the fiber optic difference. 

Needless to say that the procedure described in this subsection can be used to test 
the self consistency of the timing calibration. The use of wrong coordinates of laser 
balls or PMTs will reveal itself as mismatch in the slewing curves (t). In fact, using this 
procedure, it was discovered that coordinates of several PMTs were interchanged. An 
extension of the procedure described here is discussed in [18] where coordinates of the 
laser balls themselves are allowed to vary and can be restored. 

Correcting the slewing curves by corresponding fiber optic delays and water propa- 
gation times yields the final calibration curves of T start as a function of ToT. 
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B.3 Photo-Electron calibration. 



The main purpose of the photo-electron calibration is to find a relationship between the 
observed ToT and the number of photo-electrons emitted inside the photo-tube. The 
calibration procedure is based on a well known occupancy method described in the 
literature (see for instance [6] or [32]) and proceeds in two general steps. First, the 
ToT-PE conversion is established for low input light levels using the occupancy method 
and then a different procedure is applied to calibrate PMTs at high light levels given 
the characteristics of the calibration filter wheel. The whole procedure relies on the 
assumptions that the number of photo-electrons produced in a PMT is proportional to the 
light intensity at the PMT's photocathode, that the input light level into the calibration 
system is constant, and that all light level modulation is due to a controlled change in 
the transmittance of the filter wheel only. 

The calibration data required for photo-electron calibration is the same as for the 
timing calibration which is obtained with laser light passing through a filter at different 
transparency settings. While it is difficult to establish the light level stability of the laser 
output, it was found that if probability of the laser to produce no light when it is triggered 
is less than 2.5%, the PE calibration results are self consistent. 

This section presents the main ideas of the PE calibration followed by a description 
of innovations in the method implementation. Full description of the occupancy method 
applied for Milagro calibration is presented in [14], [32] and [13]. 

B.3.1 Low light level calibration and the Occupancy method. 

The occupancy method is based on the assumption that the number of photo-electrons 
produced at a PMT's photocathode obeys a Poisson distribution: P(n; A) = ^e~ A . 
Here A is the mean number of PEs produced at the photocathode. This is justified by the 
assumption that the emission of a photo-electron is not related to emission of a different 
one which is true if the photo-tube did not reach its saturation. 

The probability rj that a photo-tube registered the light pulse (which means at least 
one photo-electron was emitted from the photo-cathode) is called the occupancy and is 
given by: 

7] = P(n > 0; A) = 1 — P(n = 0; A) = 1 — e~ A A = -ln(l-??) 

Based on its definition, occupancy rj can be easily measured if a PMT is illuminated 
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many times by the light pulses of identical intensity: 

number of observed pulses 
number of sent pulses 

As the intensity of input light is varied, it is possible to find a relationship between 
the number of PEs and the observed ToT directly. However, for high light levels when 
A > 2, it is not possible to measure A reliably based on 77, since the error of the mea- 
surement on A increases exponentially with error on 77: 

AA = —^—An = e x An 

1-7] 

B.3.2 High light level calibration. 

The high light level calibration is based on the assumption that there is no saturation of 
the PMT channel and the mean number of photo-electrons produced is proportional to 
the input light level intensity. If the transmittance of the filter wheel T is known, then: 

A = a-T 

where a is some parameter which is constant, but different for different laser ball-PMT 
pairs. It can be found from this equation at low light levels because T is known and A 
can be measured with the occupancy method. Thus, given the transmittance properties 
of the filter wheel, the ToT to PE conversion can be found at high light levels with a 
linear error on A. 

For some PMT laser ball pairs, even the lowest light level possible was relatively 
high for the occupancy method to be used. For these PMTs, a farther away laser ball 
was used to establish the ToT-to-PE conversion for lower light levels. The obtained 
conversion curve was then extended by the data from the nearest laser ball to the highest 
possible light level. 

If, contrary to the assumption, saturation of the PMT channel is present, the number 
of PEs can not be established using this method. However, since the goal of calibration 
is to study the PMT response to different light levels, this is not a problem and the 
method described here allows one to infer the light intensity at the PMT cathode from 
the observed ToT as the effective number of PEs which should have been emitted from 
the photocathode provided the PMT response were linear. 
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B.3.3 Filter calibration. 

As was mentioned earlier, transmittance properties of the filter wheel are important for 
high light level PE calibration. The transmittance properties can be obtained from the 
manufacturer of the filter or can be measured in laboratory. A method to calibrate the 
filter wheel using the same calibration data as for slewing and PE calibration was pro- 
posed and used. This method has the advantage that the filter is calibrated as it is being 
used in detector calibration. The method employs the occupancy method with additional 
supposition that for any two sufficiently close transmittances of the filter Ti and T 2 there 
exist a PMT for which the occupancy method can be used at both light intensities. Then, 
given two corresponding measurements of mean number of photo-electrons: 

Ti _ A2 

The transmittance T 3 can be related to T 2 in the analogues manner and so forth, 
leading to the restoration of the levels of transmittance for all filter settings. Because the 
absolute calibration of the filter is not required, it is always possible to set T — 1. 

B.3.4 Dynamic Noise Suppression. 

The PE and filter calibration procedures rely on correct knowledge of PMT occupancy. 
PMT thermo-electron emission, Cherenkov light from the shower particles and other 
sources can cause a signal on the PMT output not related to the input calibration light. 
This noise will increase the measured PMT occupancy and damage the calibration ac- 
curacy. 

Dynamic noise suppression is a technique which allows correction of the apparent 
occupancy on a tube by tube basis and is based on the assumption that the arrival of 
the laser light is not correlated with the noise pulses. Then, the probability to observe 
anything (apparent occupancy P(any)) is: 

P{any) = P(laser) + P(noise) — P(laser) • P {noise) 

. P {any) — P {noise) 
n = P( laser) = — — . . 
' V ; l-P{noise) 

where P(laser) is the probability to observe laser light (true occupancy rj) and P{noise) 
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Figure B.4: Filter wheel calibration with and without noise suppression. (Bank 2,3 
represent laser balls 11-20 and 21-30 respectively, while AS and MU represent PMTs 
from "top" and "bottom" layers used for filter calibration.) 



is probability to observe noise pulse. 

P(any) can be measured by sending laser pulses to a PMT, as before, while Pinoise) 
can be estimated by sending uncorrected triggers to data acquisition system without any 
light input which can be interlaced with the laser data taking. 

The dynamic noise suppression is an important step in calibration process and was 
used for the ToT-to-PE and filter calibrations. The effect of the noise suppression is 
shown on figure B.4 where the filter wheel calibration curves are presented with and 
without the noise suppression. 



B.3.5 Statistical error of the occupancy method. 

In order to address the question of the accuracy of the occupancy method, suppose that n 
shots of the laser beam were sent, out of which m were detected by the PMT. Occupancy 
rj is estimated as: 

m 

V = — 

n 
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7] 


s(v) 


pass/fail 


0.1 


0.009 


fail 


0.2 


0.009 


pass 


0.3 


0.008 


pass 








0.8 


0.004 


pass 


0.9 


0.002 


fail 



Table B.l: Occupancy accuracy test results to satisfy the error q = ^ = 0.01 on 
measured A with n = 2000 laser shots. 

The question of the accuracy of this estimation is that of the confidence interval. 
Since the probability of a PMT registering a signal in a shot is equal to rj, then the 
number of detected pulses m is distributed according to binomial distribution: 

77 I 

P v (m,n) = OT(l " Vr m , = ■ 

m\(n — m)\ 

The upper bound of the confidence interval r] upper corresponding to the significance 
a is defined so that probability of detecting k < m pulses is less than (1 — a): 

m 

P Vupper (k < m, n) = £ C k n r) k upper (l ~ Vu PP er) n - k < 1 - OL (B.l) 

k=0 

Correspondingly, lower bound r\i ower is such that 

P Vlow Jk >m,n)=jr C k nV k ower (l - w) n " fc < 1 - a 

k=m 

or 

m— 1 

E CnVLeAl - W) n " fc > a (B.2) 

k=0 

If the required relative error on the value of occupancy is 5 {if) = ^ then r] upper 
should be no greater than (1 + S(rj))r] and r\i oweT should be no less than (1 — 5(r)))r). 
Thus, for given number of laser shots and number of registered pulses it is possible 
to check if the required accuracy is met. The requirements on the accuracy 5(rj) are 
governed by the requirement on the relative error q = on PE (A) determination: 
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A = — ln(l — 77) A\ = —^—Ar] 

1 — 1] 

A\ At] l-il Wl v 

q = — = -(i- v )Hi- v ) s( v ) = - q —Hi-v) 

Now, the task is to estimate the allowed range of occupancies 77 with error not ex- 
ceeding the specified value q provided that n laser pulses were sent. 
For 77 < 0.5 the equation B.2 will be automatically satisfied if: 

r]*n 

E CnVtpperi 1 ~ VupperT^ < 1 ~ «, Vupper = (1 + Kv))V 
k=0 

and for 77 > 0.5 the B.l is satisfied if: 

r)*n— 1 

E CnVtweri 1 ~ Viewer)^ > «, Viewer = (1 ~ <$fa))»7 
fc=0 

Therefore, for relative error g = 0.01, confidence a = 95% and n = 2000 the 
allowed range of occupancies to be used is 0.2 < r\ < 0.8 (see table B.l). For the 
phototube noise measurements (77 ~ 0.03) the requirements are eased: q = 0.1. We 
also verify that using 60000 random "shots" is just enough to reach the goal. (The 
requirement of q — 0.01 is met at 66% confidence level only.) 



B.3.6 Threshold effect on the occupancy measurement. 

The problem addressed here is that of a finite threshold which a PMT pulse should cross 
in order to be recorded by the electronics. The presence of the finite threshold leads 
to an under estimation of the occupancy and if the final effect on the number of PEs is 
bigger than q a correction should be made. 

The assumption of PMT operation is that electron multiplication in each stage is 
a statistical process with some average gain g. Then, the number of electrons k on 
the output of the first amplification stage obeys Poisson distribution with average gw 
where w is the number of photo-electrons emitted from the PMT photocathode (This 
assumption is similar to that of occupancy method and is based on proposition that 
probability of emitting an "amplified" electron does not depend on other electrons being 
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emitted.) and is distributed according to: 

It is the number of electrons on the output of the first stage that dominates the fluctu- 
ations on the output of the entire cascade and is responsible for ability of the PMT pulse 
to cross the discriminator threshold. The PMT electronic channel was constructed in 
such a way that only pg or more electrons after the first stage will produce signal strong 
enough to be detected. Thus, the probability that PMT signal of strength w PEs will not 
be detected is: 

k<P9 k<pg , s k 

/?(«,)= £p(fc;<H= E^rr^ 

k=0 k=o K - 

and the occupancy given average number of PE A is decreased: 

oo \ 2 

77(A) = 1 - P(0; A) - £ P( W ; A)/3H = 1 - e~ x - \e~ x (3(l) - -e- x (3(2) 

w=l 1 

For estimation purposes Milagro PMT is assumed to have uniform stage to stage 
amplification with typical gain of the entire PMT's cascade of 2 • 10 7 (see [4] for the 
measured PMT's gain at the operated voltage). Since the PMT consists of 10 stages, the 
gain of a single stage is g — (2 • 10 7 ) To = 5.3 and the threshold level is set at p = 0.25 
(signals with more then 0.25 equivalent PEs will cross the discriminator threshold and 
will be detected.) the function (3(w) falls off rapidly and for given PMT parameters 
a one PE input will be lost with probability (3(1) = 3.1 ■ 10~ 2 while 2PE — with 
(3(2) = 2.9 • 10~ 4 and can be neglected. Hence, taking into account only the loss of 1 
PE signals, the measured occupancy is equal to: 

77 « 1 -e~ x -(3(l)\e- x 

The magnitude of the correction (3(l)\e~ x has to be compared with the required 
accuracy on 77 which leads to the direct comparison of (3(1) and q. 

Therefore, it is concluded that the systematic error on occupancy is comparable to 
the statistical one and for the desired accuracy of PE determination (few per cent) the 
threshold effect can be neglected. If the probabilities (3(1) are known for each PMT the 
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correction, could be done with the following approximate formula: 



Also note, the filter wheel calibration is quasi-immune from this problem because it 
is based on the ratio of the A's for a given PMT and the jzjxrj f actor cancels. 

B.4 Calibration Extrapolation. 

The maximum light level at which the calibration data was available often was lower 
than could be observed in the shower data rendering strong PMT pulses unusable. To 
cope with this problem, extrapolation was used to infer the values of calibration parame- 
ters beyond the calibrated range based on the known values and trends. Indeed, typically 
it was required to extend HiToT calibrated range by about 100 (ns) to interpret shower 
data. 2 

B.4.1 Slewing extrapolation. 

The shape of the slewing correction function depends on the discriminator threshold 
level, amplification coefficients, gains of PMTs, wiring and so on. Instead of trying to 
take all the unknown parameters into account and putting together a physical model of 
the slewing, a statistical one was built taking the following approach. 

It is believed that all PMT channels (PMTs themselves and electronic boards) were 
designed and manufactured to meet common characteristics. Therefore, the study of 
the channels' responses (calibration) can be viewed as a multiple (about 700 times) 
measurement of a singe function: T start vs ToT. The fact that the curves obtained for 
different channels are slightly different can be attributed to the "manufacturing imper- 
fections" and the channels differ only due to unavoidable uncontrollable reasons such as 
spread of characteristics of electronic components and/or actual slewing measurement 
errors. Thus, a slewing curve for a PMT can be viewed as a particular realization of 
some random function. All slewing curves together form one slewing function fam- 
ily characterized by its mean dependence 3 m(t) and correlation function K(t,t'), by 

2 LoToT was never extrapolated as it is prudent to switch to the use of HiToT as soon as it becomes 
meaningful. 

3 Here, t and t' denote the time over threshold. 



91 



analogy with the mean and dispersion for a random variable. 

The characteristics m(t) and K(t, t') of the random function were deduced from 
the observed high range slewing calibration data and, based on the random function 
framework carried out to the first order of canonical expansion, the value of the slewing 
correction x(t) given the last known calibrated value x(ti) at time t x is: 

, N , \ ^(^l) _ m {t\) T r, X 

= m(t) + x( tl ,ti) X(Ml) 
with the root mean square error of: 

(#(Mi)) 2 



rms(t) 



Using this method slewing curves were extrapolated only to the point where real 
data for at least 50 PMTs existed with estimated error on extrapolation of the order of 
0.7 (ns). 4 Beyond that, linear extrapolation was used with the slope of 0.0381 jff^ - 

Details describing the extrapolation method used to extend slewing curves for the 
HiToT calibration can be found in report [16]. A brief review of the notion of random 
functions can be found in memo [20]. 



B.4.2 PE extrapolation. 

Contrary to slewing calibration where the quality of the data increases with the input 
light level, the quality of ToT-to-PE conversion degrades due to exponential relation 
between ToT and PE and possible PMT saturation. This lead to the conclusion that 
sophisticated random function extrapolation is not justified and the extrapolation was 
developed based on a simple physical argument that log PE ~ ToT. Thus, the PE vs 
ToT data was fit to a third order polynomial of the form: 

In PE = a + a{ToT + a 2 ToT 2 + a 3 ToT 3 

and the values of the polynomial were used as the extrapolation. Needless to say that 
beyond approximately PE = 100 the error grows very fast and any algorithm relying 

4 The comparison of extrapolated data from a calibration run with calibrated data obtained indepen- 
dently [38] yielded the measured extrapolation error of 0.55 (ns) in good agreement with the expectations. 
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on PE should treat the values beyond PE = 100 as logical "big", "Big" and "BIG". 



B.5 Energy calibration. 

It is planned that absolute energy calibration measurements will be done using through- 
going muons. The imaging capabilities of the detector will be exploited in order to find, 
fit and select well-defined through-going muon tracks. Once the geometry of the track is 
known, the Cherenkov energy deposit will be estimated and compared against the photo- 
electron distribution in the event. This was the primary absolute energy calibration 
method used in the 1MB detector [6]. 
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Appendix C 

Auxiliary Celestial Coordinate system 



An additional coordinate system on the Celestial sphere which is used in this work is 
defined relative to a preselected celestial object (point on the celestial sphere). The 
object is the origin of the coordinate grid. The "zero" axis of this coordinate system can 
be chosen to point to any point M on the celestial sphere. 

Let L = (Sl, oil) be the point on the Celestial sphere which should be the center of 
the map directed to point M = (5m, olm), then the point X = (5x, ctx) on the Celestial 
sphere will have coordinates (x, relative to the point L defined as (see figure C.2): 

x =LX= ILCX, f = IMLX 

£ is measured clockwise from the LM line and £ e (— n; 7r), % G (0; n) 
From spherical triangle ALPX: 

cosx = cos(LX) = siylSl sin^x + cos5l cosSx cos(a^ — osx) 

From spherical triangle APMX: 

cos(MX) = szhSm sin 5x + cos 5m cos 5x cos(q;m — osx) 

From spherical triangle APML: 

cos(ML) = sin 5m sin 5l + cos 5m cos5l cos(«m — oil) 
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p 




Figured: Schematic of the map Figure C.2: Auxiliary coordinates on the 

Celestial sphere 

From spherical triangle AMLX: 

cos(MX) = cos(ML) cos(LX) + sin(ML) sm(LX) cos(MLX) 

, , , r vx cos(AOT) - cos(ML) cos(LX) 
cos£ = cos(lMLX) = — i > - v L — ^ L 

sin(ML) sin(LX) 

There is no problem in defining the x, however, points symmetric with respect to 
MCL-plane are indistinguishable. A way to solve this problem is to introduce an aux- 
iliary vector n = cVl x Ut. 1 Then if (n ■ cl) > then £ e (0; tt) else f G (-tt; 0). 



'Since Right Ascensions of the objects come only in differences, the derived formulae will work in 
the local coordinates of (S,H) by substituting a* — > iJ*. However, the (<5, iJ) coordinate system is 
left-handed and the cross-product is defined in the right-handed system. Beware! 
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Appendix D 



Kinematics of the particles in the Solar 
system (Simulations appendix) 

D.l propagate_inf inity ( ) . 

The task of this function is to find a particle state at infinity given its state at the annihila- 
tion point by "backward-in-time" propagation. The function, however, should return the 
state in "forward-time" frame since only these quantities have physical meaning. The 
sequence of actions to be performed is summarized below: 

1 . Find orientation of the orbit in space (section D.2). 

2. Check that a particle crosses the Sun (section D.3). If "no", go to 6. 

3. Generate column density which should be accumulated before next scattering 
(section D.10). Compute parameters of the inside and outside orbits (if orbit at 
least partially exits the Sun) (see section D.5). Propagate the particle to its first 
scattering point or go to 6 if the particle does not scatter. (The last situation is 
not possible within the considered model, but is implemented for future develop- 
ment.) 

4. Generate scattering off of a proton (section D.9), compute the parameters of the 
new orbit (section D.5). Generate the column density to be accumulated until the 
next scattering (section D.10). 
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• If the particle in on unbound orbit, propagate it to the edge of the Sun, com- 
pute the Runge-Lenz vector and go to 6. 

• If orbit is completely inside the Sun find the next scattering point (sec- 
tion D.6). 

• If particle leaves the Sun, rotate the orbit according to the column density 
required and then, propagate the particle to the next scattering point (sec- 
tion D.7). 

5. Go to 4. 

6. Use Runge-Lenz vector to find the Voo or record that the particle is on a non- 
crossing Sun bound orbit. Perform time reflection on initial and final state (sec- 
tion D.7). 



D.2 Conserved quantities. 



In the considered model, the Sun is a ball of uniform density of mass M and radius R Q . 
Thus, the gravitational potential U (r) is the function of the distance r from the center of 
the Sun and is: 



U(r) 



-7 r > R G 

_a 2 _ 3a _ a ( r 2 q\ <- p 



where a = G N M & and G N is the gravitational constant. The energy and angular mo- 
mentum are conserved in such a system and are: 

v 2 

£ = — + U(r) = const, J — f x v — const 

When the particle moves outside of the Sun, its trajectory can be described by a 
conical section with the Sun in one of its foci. The trajectories of the particles inside the 
Sun are elliptical only with the center of the Sun in the center of the ellipse. 

There is an additional conserved quantity when the particles is outside the Sun. It is 
the Runge-Lenz vector IC: 



JC = fx J -ar/r |/C| = Va 2 + 2£J 2 
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The Runge-Lenz vector points along the major axis from focus to perihelion. 



D.3 Which orbits cross the Sun. 

The energy conservation law outside the Sun has the form: 

2 2r 2 r 

where v r is the radial velocity of the particle. The particle will cross the Sun if its 
minimum distance to the center of the Sun r min is less than R Q . At this point v r — 
and: 



J 2 a , V« 2 + 2£J 2 ~ a 

Train. 



o r 2 r ■ »«» oc 

A I m i n I m.m ^ 



Thus, only the trajectories for which J 2 < 2R Q (£R Q + a) will cross the Sun. 

A separate remark should be made regarding the unbound orbits which cross the 
Sun. Particles on such orbits may never cross the Sun. Indeed, if a particle is at its 
annihilation point and if it happens to be on an unbound orbit (as can be determined 
from its velocity and position vectors), the angle between its velocity and the Runge- 
Lenz vector should be acute for the particle to pass through the Sun. In other words, if 
£ > and JC ■ v < the particle will not cross the Sun. 



D.4 Rotation of a vector B around a vector L by an an- 
gle 7. 

Only the proper rotations on angle 7 are considered where the rotation is governed by 
the "right handed screw" rule around vector L, \L\ 7^ 0. The new vector B' is obtained 
from the original B by application of the rotation matrix A(-f, L). 

B' = A(>y,L)B 

The explicit from of the transformation in Cartesian coordinates is: 
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B' 



( B' x 

B 'y 
\B' Z 



a + cL 2 x 
cL x L y + bL z 
cL x L z — bL y 



cL x L y — bL z 



a 



cLyL z ~\~ bL^ 



cL x L z + bLy 
cL y L z — bL x 
a + cL 2 



( B x 

By 

V B z 



where 



cos 7 b 



sin 7 



1 — cos 7 

Z 2 



If operation of rotation of vector K, is performed around vector J, a special case 
arises when J = 0. In this situation the operation of rotation is not defined. However, 
from the physics of the situation, it follows that the rotation angle 7 is either or ir. 
Thus the rotation in this case is very simple: 



(7 == 7r) 



if\7 
else 



KJ = -JC 

it' = it 



D.5 Some facts about elliptical trajectories. 

In addition to the global coordinate system which is attached to the Sun and whose z 
axis is oriented along the Sun's motion in the Galactic disk, there are several auxiliary 
coordinate systems which are convenient to introduce. Both of the auxiliary coordinate 
systems are centered on the Sun's center and one of them is used to analyze the particle 
trajectory inside the Sun, and the other — outside. 

Capital letters for the names of the variables will represent the trajectory which is 
outside the Sun and the small letters will represent the parameters for the trajectories 
which lie inside the Sun. 

D.5.1 Equation of the ellipse. 

The equation of ellipse with the coordinate system at its center is 



a? b 2 1 — e 2 cos 2 
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Figure D.l: Elliptical trajectory inside the Sun. 



where = is the point of maximal distance from the point on ellipse to the origin (see 
figure D.l). 

The equation of ellipse with the coordinate system at one of its foci 

P 

r = 

1 + E cos $ 

where $ = is the point of minimal distance from the point on ellipse to the origin (see 
figure D.2). 

The relationship between the semimajor axis a (or A) and semiminor one b (or B) 
and the eccentricity e (or E) and the latus rectum p (or P) is: 

e 2 = 1 - b 2 /a 2 p = b 2 /a 
E 2 = l-B 2 /A 2 P = B 2 /A 

D.5.2 Ellipse inside the Sun. 

The parameters of the ellipse inside the Sun are: 
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a 



^ ( £'R P) - j£' 2 R 2 Q - aJ 2 /R & 



£' = £ + 



3a 



2R 







The trajectory of a particle inside the Sun is an ellipse since the potential energy 
varies as distance from the Sun's center squared. Moreover, the ellipse is centered on 
the Sun's center. It is always possible to choose the coordinate system in such a way 
that the OX axis is along the semimajor axis of the ellipse and OY — the semiminor 
one. Because same calculations are simpler using one parameterization of an ellipse and 
others in another, two different ellipse parameterizations are used in this work. One is by 
providing the polar angle <p measured counterclockwise from the OX and the distance 
to the point from the origin r and the other is by specifying a phase angle ip "measured 
clockwise from the OY axis" only. 1 (See figure D.l.) 

The orientation of the coordinate system is chosen by requiring that at one selected 
point (usually initial point of propagation) if the angle between f and v is acute, the 
phase angle ip of the point should be between zero and n/2 and between — it/ 2 and zero 
if the angle is obtuse (see fig D.l). 

The equations of the same ellipse in different parameterizations are: 



bcosif) 
a sin ip 



tanc/> 

r 2 = b 2 + (a 2 - b 2 ) sin 2 -ip = — 



The length s of the elliptical arc between the angles and ip 2 is: 



dx = a cos ip dip 
dy = —b sin ip dip 



(ds) 2 = (dx) 2 + (dy) 2 = a 2 



1 



a 



sin 2 ip 



(D. 



(dipf 



ds = a\Jl - e 2 sin 2 ip dip, e 2 = 1 - b 2 /a 2 



s = a [ J I — e 2 sin 2 ip dip 

hi 

Thus, the length of an elliptical arc can be expressed in terms of the elliptical inte- 



'The angle (tt/2 — ip) is called the eccentric anomaly in celestial mechanics. 
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gral: 

EE(9,k)= f Jl -k 2 sin 2 tdt, k 2 < 1 6 G [0; tt/2] 
Jo 

D.5.3 Angle between r and v. 

If current position f and velocity v are known, it is easy to find the angle (r, v) between 
the two: 

^~ ^ ^* . 

cos (r, -u) = 

rv 

If velocity £f is not known, but the the phase angle ip corresponding to the position f 
on the ellipse is known it is possible to find the angle between f and v. Inside the Sun, 
the energy conservation law is: 



- rr^ J" / J 2 

J = fxv =>- 1^71 = \f] ■ \v\ sin (r, =>- sin (r, £f) 



rv V r 2 t> 2 



^ . arcsin(^/^), ^ G [ 7t; tt/2] U[0j 7r /2] ^ tan^ > 

' ' ' ' "~ 1 7T - arcsin(^JJ), ^ e [-tt/2; 0] Ubr/2; tt] ^ tan^ < 

To find v it is enough to rotate the vector f around J by angle (f, v) (see figure D.l) 
and rescale appropriately: 

v = -rotate_vector(r, J , (r, v)) 
r 

D.5.4 Ellipse outside the Sun. 

The parameters of the ellipse outside the Sun are: 

E 2 = l + ^f- P = J 2 /a 
cr 
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Figure D.2: Elliptical trajectory outside the Sun. 



The bound trajectory while outside the Sun is an ellipse with the center of the Sun 
at one of its foci ("the outside ellipse"or "the outside orbit"). It is interesting to find the 
angle $ between the Runge-Lenz vector (the direction on the perihelion) and the point 
where the outside ellipse intersects the Sun (see figure D.2). This angle is easily found 
from the equation of the outside ellipse: 

P 

— = 1 + E cos $ 

Rq 

cos $ = 1 (P/R & - 1) = J2 ~ l R& (D.4) 



D.5.5 Rotation of outside orbit due to passage through the Sun. 

The angle of rotation of the orbit in a single pass through the Sun is defined as the angle 
of rotation of Runge-Lenz vector due to this passage. 

If 7 is angle of rotation of the orbit as it passes through the Sun, and $ is the 
coordinate describing the point of entrance of the orbit into the Sun and O is direction 
to the same point, but with respect to the orbit which is inside the Sun, we obtain the 
relationship for the angle of rotation of the orbit (see figure D.3): 

0o = $o + (tt/2 - 7/2) =>• 7 = 7t + 2$ - 20o 
From the equation (D.4) $ can be obtained: 
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Figure D.3: Rotation of the outside orbit due to passage through the Sun. 



J 2 -aR & \ 
$0 = arCCOS ( _ ^o~J 

From the equations (D. 1) and (D.3) the phase and spatial angles of the entrance point 
of the orbit into the Sun are: 




. . - b 2 ftcos^o 

sniV>o = -y— tt tan0 o = — : — - (D.5) 

tr a sin ip 

The "— " sign is selected since at the entrance point into the Sun (r • v) < always. 
Rotation of the Runge-Lenz vector )C happens on the angle (—7) around J in one 
passage of the orbit through the Sun. 



D.6 Motion inside the Sun propagate_in_sun ( ) . 

Only the motion which is initialized from inside and while inside the Sun is treated here. 
That is, the task considered here is knowing the requested column density to evolve a 
particle from its initial position and velocity inside the Sun to the scattering point inside 
the Sun or if the particle exits the Sun — output the new direction of Runge-Lenz vector 
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and the remaining column density to be accumulated in subsequent passes through the 
Sun. 

There are several situations possible: 

• The orbit is completely inside the Sun. Propagate the particle to its next scattering. 
Output final f and v. 

• The orbit exits the Sun and column density to be accumulated is large to allow the 
particle to exit the Sun. Output the final Runge-Lenz vector )C. 

• The orbit exits the Sun geometrically, but the column density to be accumulated 
is not large enough. Particle scatters before its exist. Output f and v just before 
the scattering. 



D.6.1 Particle is inside the Sun. 

An orbit lies completely inside the Sun if the semimajor axis of the ellipse is not larger 
than the radius of the Sun: a < R G . The phase angle ipi and the spatial angle 0i 
corresponding to the current position fi inside the Sun can be found from the equation 
of the ellipse in phase coordinates (see equations (D.3) and (D.l)): 



sin^i = sign(fi ■ v^J V'i ^ h 71 "/ 2 ; V 2 ] 

(D.6) 

tan & = ^1 

The final angles ip 2 an d 02 are such that the column density accumulated by a par- 
ticle traveling inside the Sun is equal to the specified column density L. Currently, the 
equation formulated in the phase angles is considered: 

L = tra jectory_length('0 1 , ip 2 ) 
This equation can be solved (see section D.8) and the corresponding spatial angle is: 

b COS 1p2 



tan 02 



a sin -02 



Now, the sought for direction of f 2 can be found by rotating the vector fi on angle 
(02 — 0i ) around J and the magnitude of r 2 an be found from the equation of the 
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ellipse (D.l): 



r\ = b 2 + (a 2 — b 2 ) sin 2 i[) 2 r 2 = rotate_vector(fi, J ', <pi — 4> 2 ) • 




The velocity vector v 2 at the point f 2 can be found by rotating the vector f 2 on the 
angle between f 2 and v 2 obtained from the equation (D.3). The magnitude of v 2 is found 
from the energy conservation law (D.2). 

In the case of a circular orbit a = b and the previous logic will fail because the initial 
phase ipi is arbitrary and the most appealing choice is to set ipi = 0. The logic of the 
presented algorithm is intact. 

D.6.2 Orbit exits the Sun, but the particle doesn't. 

It is possible to have a situation where a particle has its initial position inside the Sun, 
but its orbit leaves the Sun geometrically. It is also possible, however, that the column 
density to be accumulated is small and the particle scatters before it has a chance to 
exit the Sun. This is case of motion confined to the interior of the Sun and was treated 
above. To make sure that the particle stays inside the Sun, the column density which can 
be accumulated until the exit from the Sun should be greater than the requested column 
density L: 



where ipi is the initial phase defined in equation (D.6) and is the phase of the exit 
point from the Sun and computed as in equation (D.l): 



If it is found that the requested column density is greater than maximum possible in 
the configuration, the particle leaves the Sun and this case is treated right below. 

D.6.3 Particle exits the Sun. 

Again, as before, <pi, tpi describe the initial point (equation (D.6)) and 02,^2 is the final 
point of propagation — the point of orbit exit from the Sun (see equations (D.l) and 



L < tra jectory_length('0 1 , ipR Q ) 




(D.7) 
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(D.7), iP 2 = ^rJ: 



. , Rr.-, — b 2 , b cos ib 2 ^ ON 

sin^ 2 = +i/ ° tan0 2 = , ] 2 (D.8) 

V a 2 — (r a sin ^ 2 

The accumulated column density inside the Sun between angles ipi and ip 2 should 
be computed and subtracted from the remaining column density. 

The spatial angle between the initial point f\ and the Runge-Lenz vector /C consists 
of the angle between the initial point and the exit point from the Sun (0 X — <p 2 ) and 
the angle <3> 2 between the exit point from the Sun and the Runge-Lenz vector which is 
found from equation (D.4). Therefore, the vector r x should be rotated by (0i — <p 2 — $ 2 ) 
around the J with rotate_vector(r 1 , J", (0 X — <p 2 — $ 2 )) to find the direction of the 
Runge-Lenz vector. 



D.7 Motion outside the Sun 

propagate_outside_sun ( ) . 

The task of this function is given the direction of the Runge-Lenz vector and the column 
density column_density to be accumulated inside the Sun until the next scattering 
point find the particle position and velocity at the scattering point. Needless to say, 
the function should process only the particles which pass through the Sun. The other 
parameters which are expected to be available are the parameters of the orbit inside the 
Sun and the phase angle ip 2 defined in the equation (D.8). 
There could be several distinct cases: 

1 . The trajectory is an ellipse 

2. The trajectory is a hyperbola or parabola (but the particle crosses the Sun) 



D.7.1 Find the scattering point. 

The parameters of the trajectory do not change between scatterings and only rotation 
of the whole orbit is possible as the orbit passes through the Sun. The angle 7 of orbit 
rotation due to a single pass through the Sun was obtained in section D.5.5. The number 
of required passes through the Sun can be be determined from the required column 
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density to be accumulated and the column density accumulated in a single pass through 
the Sun. This will give the total angle of the orbit rotation before the scattering. 

Due to spherical symmetry of the Sun, the column density accumulated in one pass 
through the Sun has the form of (see section D.5.2 and equation (D.7)): 



L 1 — 2 • tra jectory_length(0, 4>r q ) 



The total angle n ■ 7 of rotation of the orbit is: 

n ■ 7 = f loor (column_density/Li) -7 



column_density := column_density — nL\ 



From figures D.2 and D.3, one sees that the vector pointing to the entrance point into 
the Sun can be found by rotating the Runge-Lenz vector by angle (— $ ), that is why 
the total rotation of IC to be performed is on angle (— 717 — $ ) around J. 

The remaining problem is to find the angle between the point of entrance into the 
Sun and the scattering point using the remaining column density to be accumulated. The 
particle will spend remaining time on an ellipse inside the Sun, and the phase angle at 
the scattering point fa can be found given the remaining column density and the phase 
angle fa (equation (D.5)) of particle entrance into the Sun (see section D.8): 

b cos fa 

column_density = tra jectory_length(wo, fa) =>• tan0i = — 

asm fa 

— * 

Now, the sought for scattering point fi can be found by rotating the vector JC on 
angle (—727 — $ + fa — fa) around J and the magnitude of fi an be found from the 
equation of the ellipse (D.l): 



r\ = b 2 + (a 2 — 6 2 ) sin 2 fa T\ = rotate_vector(/C, jf, — ny — $ O + O — fa) • 



Ni«i 



The velocity vector v\ at the point r x can be found by rotating the vector fi on 

— * 

the angle between fi and v\ obtained from the equation (D.3) around vector J. The 
magnitude of Vi is found from the energy conservation law (D.2). 
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Figure D.4: Finalize 



D.7.2 Particle on an unbound orbit crossing the Sun. 

This situation can occur if a particle is generated on a hyperbolic or parabolic orbit, 
but it passes through the Sun. If the requested column density is smaller than can be 
accumulated in one passage through the Sun, the above logic will suffice. If however, 
the requested column density can not be accumulated in a single pass, the particle will 
leave the Solar system, thus, in the previous calculations the angle of orbit rotation due 
to passage through the Sun is equal to 7 from the section D.5.5 only. The output of 
the function in this case is the final direction of the Runge-Lenz vector and is found by 
rotating the input Runge-Lenz vector on the angle (—7) around J. 

D.7.3 Velocity at infinity. 

This is the final function and its task is to find the velocity of the particle at infinity. 
Since only "forward" in time information is needed, this function should perform time- 
reflection on the initial and final velocities of the particle. 

The time reflection can be performed in the spherical coordinates by requiring the 
change in the polar angle 9 — > (n — 9) and in the azimuth — > (it + 0)modulo(27r). 
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A particle can exit the Solar system only when it is on a hyperbolic or a parabolic 
trajectory when E > 0. The direction of the velocity at infinity can coincides with 
the direction of the trajectory asymptote. Thus, the velocity direction can be found 
by rotating the Runge-Lenz vector to point to the asymptote (see figure D.4) by angle 

(it — arctani?/A). 

B J 28 



A y^2S a a v 
and the magnitude of the velocity can be found from the energy conservation law: 

Voo = V2S 

The above works for the Solar system escape on a parabolic orbit too, because in 
this case the particle starts its fall into the Solar system in (/C) direction. Even though 
Voo = in this case, it is necessary to know the direction in which the particle started 
to fall onto the Solar System, that is why should be stored in spherical coordinates 
providing the information about the magnitude and the direction of the velocity even if 

| ^oo | = 0. 

If E < and the particle does not cross the Sun, it stays on closed bound orbit 
around the Sun and the velocity at infinity is not defined. 



D.8 SolvePath4Psi. 

The task of the function is given initial point on an ellipse inside the Sun via its phase 
angle ip in find the point ip out such that the ellipse pathlength from ip in to ip out is equal 
to the given value of path_in. If such a point is found, its phase ip out is returned. 
If the point can not be found because the particle leaves the Sun before the requested 
pathlength is accumulated, the path_in is decreased by the path traveled in the Sun 
and ipout is set to the exit point from the Sun. 

The length of an ellipse from point ip = to ip is found by the elliptical integral: 

EE(ip,e) = sign(ip) ■ f \Jl — e 2 sin 2 tdt 

J 

where e is eccentricity of the ellipse and < n/2. 

Thus, elliptical distance between two points on ellipse is 
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L(ipim ipout) = a(EE(ip out , e) - EE(ip in , e)) 



The problem with this definition is that both initial ip in and final ^ out phases should 
be less than n/2, otherwise, the definition of the distance on the ellipse must be modified. 
The initial phase ip in is within allowed range by its construction. 

D.8.1 Bracketing the root. 

If the orbit exits the Sun, both ip in and ip out are within allowed bounds since the phase 
for the exit point ip out < tp R < n/2. 

If the pathlength path_in is larger than the path from ip in to ip R the particle exits 
the Sun. The path L(ip in , ip R ) is subtracted from path_in and the remaining path_in 
and ip R are returned with a flag that the particle exited the Sun. 

Otherwise, the particle stays inside the Sun even though, its orbit geometrically exits 
the Sun and ip in < ip out < tp R . The root bracket is found. If the ellipse is completely 
inside the Sun, more things have to be done. 

If the total requested pathlength path_in is greater than the length of ellipse, parti- 
cle will make one or several full revolutions which are not interesting for us. Pathlength 
path_in should be reduced by the path acquired in full revolutions. After this is done, 
ip out is within: ip in < ip out < 2n + ip in . The problem is reduced to the one where the 
particle does not make a full revolution in the Sun. 

Since it is known that the particle makes less than a full revolution, a check can be 
made if %p out < n/2 and ip out < 3n/2 by comparing path_in with L(ip in ,n/2) and 



If it is found that the solution is in ip out < n/2 the lower bound on the root ipi is set 
to ipi n and the upper ip u to n/2. 

If it is found that the solution is in n/2 < ip out < 3n/2, the pathlength L(ip in , n/2) is 
subtracted from path_in and ip in is reset to to —n/2, ip t is set to —n/2, and the upper 
bound ip u is set to n/2. This effectively produces a rotation of the coordinate system by 
the angle of n. Thus, after the root is found, its value will have to be increased by n. 

If it is found that the solution is in %p out > 3n/2, the pathlength L(ip in) 3n/2) = 



in , n/2)+L(—n/2, n/2)) is subtracted from path_in and becomes the upper 



bound on the root ip u = ip in . The lower bound on the root tpi is set to —n/2 and ip in is 



(L(^ n ,7r/2)+L(-7r/2,7r/2)). 
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reset to to —ir/2. This effectively produces a rotation of the coordinate system by the 
angle of 2ir. Thus, after the root is found, its value will have to be increased by 2ir. 
When the algorithm arrives to this point, the bracketing of the root is finished: 

— 7r/2 <ipi< ipout < ipu < 7r /2 with possible flag to increase ?/w by ir or 2n. 

D.8.2 Circular bracketing. 

The phase angle interval can be reduced further by noticing that the length of the arc of 
radius of semimajor axis a is not smaller than the length of the ellipse within the same 
boundaries on the phase angle ip. Likewise, the length of the arc of the semiminor axis 
radius b is not greater than the length of the ellipse within the same boundaries on the 
phase angle ip. In other words, the solution ip out can be bracketed as: 

^i<^out<^2, ipi = (ipi + 2n/a) ip 2 =rran((^i + 27r/&),^ u ) 

D.8.3 The solution. 

The equation for t/w which needs to be solved is: 

L{ipin, ipout) = path_in 

This equation can be solved on the specified interval (ipi, ip 2 ) using the secant or 
bisection method until the error on accumulated pathlength L becomes within allowed 
range. (It was found that for e > 0.99, when the ellipse is close to degeneration into a 
line, the bisection method is faster than the secant.) 

The solution of this equation should be increased by n or 2n if coordinate system 
rotation was required as described before. 

D.9 Scattering in the Sun. 

Let momenta of a neutralino and a scatterer (proton) before scattering in the laboratory 
reference frame be p x and p p and after p' and p' p accordingly. It should be noted that 
two vectors p x and p' form a plane thus p x can be found from p ' by a rotation. The 
magnitude and the angle of rotation can be found from the kinematics of the elastic 
scattering. The axis of rotation has a random direction in space. 
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Figure D.5: Scattering diagram. AO 

— * 

OC = \fxv\n. 



m x +m p 



(Px + Pp), OB = ^p^(p x + P P ), 



D.9.1 Elastic scattering. 

In the center-of-mass frame, the momenta before scattering are: 



P x o = V v 
PpO = -fJ>v 



m x m p 
m x + m p 



Px 



m x m p 



r] + l 



\Px ~ VPp\ V = 



m r , 



In the center-of-mass reference frame, the act of elastic scattering can change the 
directions of the momenta only 2 : 

Pxo = l*vfi = ^[\p x ~ VP P \n 

I Ppo = -pvn = -^TilPx - VP P \n 

where n is a unit vector along the velocity of the neutralino after the collision. In labo- 
ratory reference frame: 



' p' x = ^vn + ^^(p x + p p ) = ^\p x - V p p \n + 



r?+l 



m x +m p 



(Px + Pp) 



P 



-^n + ^^- p (Px + Pp) = -^i\Px - VP P \n + 



m x +m v 



(Px + Pp) 



2 This is merely a statement of the energy conservation law 
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The task is to find the p x when p ' is known. If in the laboratory frame the scatter is 

at rest before the scattering p p = 0, then (see figure D.5): OB = /iv AB = p x and point 
B is on the sphere. In this case, the angle of deflection of neutralino velocity 9 1 from its 
original direction is related to the scattering angle 9 in the center-of-mass frame as: 

m v sin 9 sin 9 

tan 0i - p 



m x + m p cos 9 rj + cos 9 
and 



IPxl = \p x 



m x + m p _ _ rj + 1 



y — \y _ _ 

+ m 2 p + 2m x m p cos 9 X x Vv' 2 + 2r]cos9 + l 



12 _ i.-?/ i2 ^ + 2r?+l 



\v x \ = \v 



x rf + 2rj cos + 1 



The energy loss in a collision is: 



(Vxl 2 -I^xl 2 ) _ 2r ? (l-cos0) \v x \ 2 
2 (1 + r]) 2 ' 2 



D.9.2 Choosing the axis of rotation. 

Vectors v x and v' x from a plane which is defined by some vector n. A vector n perpen- 
dicular to v' x can be found by solving: 

v' x • n — v' x n x + v' y ny + v' z n z = 

This equation can be solved by setting trial coordinates for the vector n — (1,1,1) 
and then modifying one of the coordinates to satisfy the orthogonality condition: 



if (u^O) n z :=-(v' x + v' y )/v' z 

else if (v' y ^ 0) n y := -v'Jv' y 
else n x := 
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Then, the obtained vector n can be rotated on a random angle around v' which will 
produce a vector perpendicular to v' and pointing in a random direction in space. The 
obtained vector defines the (v ' , v x ) plane. 



D.10 Generate path inside the Sun. 

The probability dp that a particle will travel distance x through matter without scattering 
and then scatter immediately after that in the distance (x, x + dx) is: 

dp = \e~ x/x dx 
A 

where A is the mean free path of the particle in matter. 

Thus, the pathlength which need to be accumulated in the Sun until next scattering 
should be drawn from an exponential distribution with parameter A. 



n p is the concentration of the scatterers and a vx is the crossection of the scattering 
process. In the given model n p = i^i m , ™ v is the mass of scatterers (protons). 
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Appendix E 



Comments on the upper limit 
construction procedure 

Very often physicists consider the question of validity of a new theory to be equivalent 
to a non-zero value of the parameter(s) of the theory. Thus, often, the tests of validity of 
a new theory are designed in such a way as to measure the value(s) of its parameter(s). 
If the measured value is non-zero it is concluded that the new theory is valid with the 
obtained value of the parameter. Such a "physical" approach seems to give adequate 
results but is not correct from a methodological point of view. Indeed, one can always 
assume that some theory is correct and obtain a non-zero value for the parameter of 
the theory based on the experiment, but that does not mean that the theory correctly 
describes the observed process. Also, if the measurement was "not successful", that 
is, the experiment could not show that the value of the parameter is non-zero, it is not 
possible to decide if the new theory is valid or not. And more importantly, the new 
theory might not even have a free parameter to be measured. 

In the defining work by Neyman [39] on statistical estimation the question of a sta- 
tistical test is separated from the question of the measurement. The procedure for pa- 
rameter estimation (measurement) demands that it is known that the process, whose 
parameter is being measured, exists and the observed data is described by the known 
distribution with the parameter being measured. If it is not the case, the procedure can 
not guarantee that the constructed confidence interval will contain the true value of the 
parameter with requested probability. This fact if often overlooked. 

When a new theory is proposed the experiment should not try to estimate the value of 
the parameter of the theory, but, instead, it should be designed to exploit the differences 
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between the adopted (old) and the new theories and check if indeed there is evidence 
to reject the old one. It is the difference between the old and the new theories which 
provides the "signal" in the test. The test of the validity of the new theory should be 
designed in the spirit of the proof by contradiction (also called indirect proof) method. 
In other words, it should be assumed that the new theory is wrong and the old one is 
correct. If the contradiction between the old theory and observed data is found (i.e. the 
"signal" is found), the assumption of validity of the old theory should be rejected and the 
new theory may be accepted as a valid one. If the contradiction is not found, nothing can 
be said regarding the validity of the new theory and the old theory can not be rejected in 
favor of the new one due to lack of evidence. 



E.l Sensitivity and upper limit. 

There is a question which an experimenter should answer when designing the experi- 
ment: what is the probability to accept H due to pure chance when H\ is true i.e. what 
is the power of the constructed test. The power of the test depends on the alternative 
hypothesis and its parameters as well as on the null hypothesis. If a new theory provides 
a large power of the test and is true, the contradiction between the observed data and the 
old theory will be found easily by the constructed test. If, on the other hand, the new 
theory provides small power if it is true, the contradiction between the observed data 
and the old theory will not be found easily. It is seen that the "strength" of a "signal" is 
defined in terms of the power of the test. The new theory predicts a "strong" detectable 
"signal" if it provides large power of the test. 

Suppose that the alternative hypothesis has the form of p±(x; A) where x is observed 
quantity and A is the parameter of the new theory. For some values of A the power 
of the test will be small and for others the power will be big. It is proposed to define 
"sensitivity" of the test (or experiment) as such values of the parameter A for which 
the test has the power of at least 50%. The sensitivity defined this way has several 
important properties: it is a detector feature and can be estimated before the experiment 
is performed, it is not a random number and does not depend on the value of the observed 
quantity. 

It is also proposed to state the upper limit when H is not rejected as such values of A 
for which the power of the test is big (say at least 90%). The choice is motivated by the 
logic that it is hard to miss such a strong "signal" and yet it was not observed. In other 
words the values of A for which there is a big chance to not reject the null hypothesis are 
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below the upper limit and the values of A for which the chance to reject H is big can be 
dismissed when H was not rejected based on the observed data. 



E.2 Problem with the current approach. 

Currently, the upper limit on the value of the parameter A of a theory is reported as 
the upper boundary of the one-sided confidence interval on the parameter which is ob- 
tained assuming that the data came from a distribution characterized by pi(x; A). This 
is, however, true and correct only in the cases when the physical process originating 
the observed data exists and according to the theory is described by pi(x; A). Again, 
only if there is no question of existence of the process is it correct to use the procedure 
for confidence interval construction for parameter estimation. Note that the bounds of 
the confidence interval are random by nature because they are functions of a random 
variable and will vary with the data observed. 

The application of the same technique to construction of a confidence interval on 
A due to a process which is not known to exist when the hypothesis of absence of this 
process is not rejected based on the observed data will lead to meaningless results. The 
constructed confidence interval will not have the desired confidence level and more over, 
the theory H which has no parameter A has not been rejected based on the observed 
data. It is irrational and illogical to assume the validity of Hi when H is not rejected. 



E.3 Example I. 

Let us consider a situation when according to the old adopted theory a quantity X is 
distributed according to Gaussian law with zero mean and known standard deviation a: 

If a new theory is correct, the data X should be distributed according to the Gaussian 
distribution with positive mean A: 

Pi(x;m) = -^ =e -(^) 2 /^ 2 A > 
Let us assume that if x > 3a the old theory is rejected and a confidence interval on 
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the values of A may be constructed. If X < 3a the old theory is not rejected and an 
upper limit on the values of A should be set with the confidence of 90% (or 1.28er). 

Further, let us assume that the outcome of the measurement is x = 0, which is the 
most probable outcome when the old theory is true. According to the current approach, 
the upper limit on the A would be A < 1.28a (This is a standard one-sided confidence 
interval on A with 90% confidence.) 

Now, suppose that the new theory is true with the parameter A = 1.29a. This 
signal is above the upper limit. What is the probability of discovering the signal in this 
experiment or what is the power of the test? For a discovery to happen x > 3a should be 
observed which has a probability of happening of 5% only. Therefore, the experiment 
which set a limit of A < 1.28<r has almost no capability to discover a stronger signal. 
This is clearly unsatisfactory result. 

Using the approach proposed here, it is required to state the upper limit correspond- 
ing to 90% power of the test. Since in the observed data x < 3a. no discovery is made. 
The upper limit would be A < (3.0 + 1.28)er = 4.28cr. The result is stated as: the upper 
limit corresponding to significance 1.35 ■ 10~ 3 and the power 0.9. 



E.4 Example II. 

The importance of knowledge of existence of a process before a measurement can be 
performed can be illustrated on a somewhat artificial example. Suppose according to 
a new theory elephants have wings. An attempt to measure the length of the wings, 
for example, resulted in the limit that their length is between and 2 centimeters, for 
example. The approach proposed here would state that no wings were found, but if they 
were longer than 5 centimeters, for example, they would have definitely been found. 

This illustrates that by following the assumption that the new theory is correct with- 
out testing it, one is in danger of reporting a limit on an absurd parameter while ac- 
cording to the approach proposed it will be clearly stated that the new effect was not 
discovered and how strong the effect should have been to be discovered. 
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E.5 Example III. 



Suppose it is known that observed data X comes from a Gaussian distribution with 
unknown mean A and known variance a: 

1 (rc-A) 2 

p(i;A) = 75P e ^ 

If it is desired to estimate the value of A by a 90% one-sided confidence interval, 
given that the observed data x is x = 0, the interval on the values of A is: — oo < A < 
1.28(7. 

Suppose, later, an existence of a new process is proposed according to which the 
same data X on the same experiment should have the distribution of: 

1 (^A-M) 2 

M^; A, /i) = -^==e ^ , fi>0 

One is tempted (and according to the current approach this would happen) to state 
that based on the same data the total signal is — oo < (A + ji) < 1.28a. This is, of 
course, incorrect, because it is not known if the new process with the parameter ji exists. 
Instead, one should consider a statistical test with: 

Hi. Pl (x; X,fj) = 7=fe" ( 2 ^ M) , /i>0 

If the value of A is known, the test would be very simple: if the observed data is 
greater than some x the null hypothesis should be rejected. For a 1.35- 1CT 3 significance 
the critical region is defined as: 

x > x = A + 3(T 

If the observed data is inside of the critical region, a confidence interval on the value 
of pi may be constructed. If the null hypothesis is not rejected based on the observed data, 
the upper limit fi u corresponding to the power of 0.9 can be found from the equation: 

/ _^ =e -{^-^?/^ dx= _^ e -vV^ dy = 0.9 

J3a+\ V2"7r<7 2 ha-iiu V27RX 2 
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If all the information available on the value of A is in the form — oo < A < X u with 
100% confidence, one is forced to construct a conservative test as above assuming that 
A = X u . This will lead to the critical region corresponding to the significance 1.35 ■ 10~ 3 
constructed as: 

X > X = \ u + 3(7 

Again, as before, if the observed data is inside the critical region, the null hypoth- 
esis of absence of the secondary process fi can be rejected and the value of fi may be 
estimated. If the observed data is outside of the critical region, the null hypothesis is 
not rejected and an upper limit on the value of \x based on power of the test can be con- 
structed as above. Note that the information on the value of \ u should come from an 
independent experiment otherwise the constructed test and upper limit are not correct. 

If, however, no information on the value of A is available the proposed test will not 
yield any meaningful result as it is not possible to distinguish data X originating due to 
process A or /i. 

A more general statement can be made. If the null and the alternative hypotheses on 
the origin of observed data X come from the same family of probability distributions 
i.e: 

Po(x; A) = f(x; A) and p ± (x; A, /i) = f(x; \ + fJ,) 

and if no knowledge regarding the value of A is available, it is not possible to construct 
a meaningful statistical test to differentiate the two. 1 However, if several independent 

'Indeed, the hypothesis test should be constructed in such a way so that the critical region on the 
values of X does not depend on A. A procedure for constructing such a region was proposed by Neyman 
and Pearson [40] for the case when: 

rflnp (x;A) d<t>(x; A) R ,, u , 
^ ^ = dX — dX — = ( ' B{X)(j}{x\ X) 

For these conditions to be satisfied the function po(x; A) should come from exponential family of the 
form: 

f(x;X) = e z W T (x)+QW+s(x) 

where Z(X), Q(X), T(x), S(x) are some functions and Z(X), Q(X) are infinitely differentiable functions 
of A. The equation of hypersurface cp(x; A) = const in the AC-space is equivalent to the equation T(x) = 
const. 

Since in the considered case both p and pi are of the same type, the likelihood ratio pi/po is inde- 
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observations x-± and x 2 can be made such that it is known that data x 1 originated due to 
process with parameter A only and x 2 is due to possible (A + /i) the test constructed in 
the form: 

p (xi,x 2 ; A) = f(xi, X)f(x 2 ; A) and p i (x 1 , x 2 ; A, /i) = f(xi; X)f(x 2 ; A + /i) 

will have the power to differentiate whether or not the data x 2 came from a new process 
even if no information on A is available before the test. 



E.6 When is the new theory valid? 

Any theory of a physical process should be considered "admissible" if no observed 
data contradicts predictions based on the theory. However, in statistical tests it is only 
possible to reject a theory in favor of some other theory and it is not possible to state with 
absolute certainty whether or not the given theory describes the given process correctly. 

The question of how much evidence contradicting to validity of the old theory in fa- 
vor of the new one should be observed in order to declare that the new theory correctly 
describes the given physical process is a philosophical one. Obviously, the answer to 
this question is of great practical importance and will govern the choice of desired sig- 
nificance and power in formulating and performing the test. 



pendent of the values of x on the hypersurface <j)(x; A) = const and the equation for the critical region 
corresponding to the error of the first kind £ becomes: 

' £ h{x;\)=constM^ = J 0(;E;A)=const: Pl/po>q Po(x; X)dx 

(E.l) 

Pi(s;A,m) _ const-{Z(\+ l i)-Z(\))+Q(\+^)-Q(\) > 

Since the integration in (E. 1) is performed over all x for which p\ /p > q is true and because pi/po > q 
is either always true or always false for all values of X, the equation (E.l) becomes: 

£ / p (x;X)dx= / p Q (x;\)dx 

•J (p(x\\)— const J (/)(x;X)— const 

For a given £ this equation can be satisfied only if po{x; A) = or if £ = 1. The former condition is 
not interesting and the later one states that the null hypothesis should be rejected all the time regardless 
of the value x. 
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